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ABSTRACT 


Cyclic steam stimulation (CSS) is the commercial oil recovery process of choice for the 
Cold Lake oil sands. Much of CSS is performed under fracture pressures and, as a result, 
performance predictions are difficult. This work presents a new mathematical model, 
involving analytical techniques, for CSS performance prediction. The model is based upon 
a fracture heating computation, coupled with fluid flow—both during steam injection and 
oil and water production. Thus two-phase flow is accounted for. Two situations, involving 
different flow geometries, bracketing actual flow in the field are considered. These consist 
of fracture flow, and flow in an elliptical geometry with a circular well in the centre. 


A vertical fracture is assumed to be created during steam injection. Fracture half length 
is calculated either from an analytical formula or provided by the user. No mechanistic 
approach is taken to do any stress analysis. An investigation for a completely closed 
fracture during production is performed. Based on the results from this investigation, 
existence of a partially open fracture during production is suggested. 


An improved average temperature calculation technique and a new method of evaluating 
the effect of spatial oil viscosity distribution have been incorporated into the new model. 


The model is used to simulate previously reported results showing an excellent match. 
A sensitivity study for the variation of different parameters produces a reliable response 
from the model. It is capable of predicting various production performance parameters ( 1.e. 
calendar day oil rate, oil-steam ratio, water-oil ratio, oil recovery etc.) of the CSS process. 
The water-oil ratio prediction is restricted to the first five cycles only. The model can serve 


as a valuable adjunct to numerical simulation or physical modelling. 
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Figure 6.74 Cumulative water produced for different steam injection 
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NOMENCLATURE 


major axis of ellipse (m) 

area (m2) 

width of fracture (m), minor axis of ellipse (m) 
specific heat (kJ m-3K-)), fluid loss coefficient (m d-!/2) 


compressibility (kPa-!), half length of the line source (m) 


diffusivity length (m) 

energy stored in the cap and base rock (kJ) 
acceleration due to gravity (9.81 ms -2) 
formation thickness (m), specific enthalpy (kJ kg-!) 
latent heat of vaporization (kJ kg-!) 
dimensionless heat input 

steam injection rate (m3 day -!) 

absolute permeability (m2) 

thermal conductivity (kJ m=! d-! K-!) 

hydraulic conductivity (ms ~!) 
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mobility ratio, volumetric heat capacity (kJ kg-!) 
pressure (kPa) 

pressure at time zero for the exponential decline curve 
fluid flow rate (m3 day ~!) 

heat loss rate (kJ s -!) 

heat flow rate (kJ day ~!) 

radius (m) 

radius in polar coordinates 

irreducible water saturation (fraction) 

oil saturation (fraction) 

residual oil saturation (fraction) 

spurt loss quotient 

initial water saturation (fraction) 


residual water saturation (fraction) 
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time (day) 
temperature (°C) 
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yu" 


U unit step function 
v variable of integration 
V volume (m3) 
Ve __ fluid leak-off velocity (m3 s-!) 
Xj steam quality 


Yo distance between the line source and the wellbore 


Greek Symbols 


thermal diffusivity (m2 s-1), effective fracture width (m) 
1, thermal expansion coefficient ( K-!) 

mobility (m2 Paz! s-!) 

fluid heat loss due to production (fraction) 
eccentricity of the ellipse (fraction) 

potential 

porosity (fraction) 

density (kg m-3) 

dimensionless x co-ordinate 

dimensionless z co-ordinate, hydraulic diffusivity 
dynamic viscosity (Pa.s) 

dimensionless y co-ordinate 

density (kg m7) 

dimensionless time 


Seco gtere oo rast OG, to), 


Ts dimensionless time for steam advancement in the fracture 
tsa dimensionless time for the steam to move into the formation 
0 dimensionless temperature, angle 


Vv kinematic viscosity (m2 s-!) 
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subscript 


0 interface between reservoir and the adjacent strata 
1 hot 

2 warm 

b bulk, interface 
Cc condensate 

d drainage 

e external 

f formation, fluid 
h hot, interface 

i initial, injected 
e) oil 


om mobile oil 


p pore 

r radial 

R reservoir 

S steam 

SC standard condition 

st stock tank, steam 

Ww water, wellbore, warm 


wc connate water 
wf flowing wellbore 


wm mobile water 


Z vertical 

Abbreviations 

CDOR calendar day oil rate 

CSS cyclic steam stimulation 

CWE cold water equivalent 

EFW elliptic flow into a well 

LGTB linear gravity with total balance 
LGWB linear gravity with water balance 


LSL large-small-large 


OSR oil-steam ratio 
WOR water-oil ratio 


Chapter I 


INTRODUCTION 


For reservoirs with medium to high viscosity oil, thermal methods of enhanced oil 
recovery are very effective. Steam was used for oil recovery as early as 1931 when it was 
injected into a sand at a depth of 380 ft. for an extended period in Woodson, Texas. After a 
relatively quiet period of 20 years with some smaller projects starting in California, the first 
large scale steam drive was employed in Tia Juana, Venezuela (Prats, 1986). Cyclic Steam 
Stimulation (CSS) or steam soak was discovered rather accidentally in the Mene Grande 
Tar Sands, Venezuela in 1959 (Adams and Khan, 1969). When steam broke through the 
overburden after injection into a non-productive primary well, the well was produced to 
relieve the pressure. To everyone's surprise, it produced a considerable amount of oil with 
little water. Since then, steam soak has been used widely all over the world. It is being 
used successfully in California where many reservoirs contain high viscosity crude oil. 
Success with this method has been reported also in Cold Lake, Alberta and Venezuela. For 
economy and high early cycle production rates, it has become routine to employ CSS for 
moderate to heavy oil reservoirs prior to a steamflood. 


The extensive use of CSS by the oil industry points to the need for a simple 
performance prediction technique. The CSS process has complex mechanisms: the reversal 
of flow during the injection and production periods: saturation hysteresis; possible 
formation of fractures during injection and, in some cases, dilation and compaction. These 
factors make the solution of simulator equations difficult, expensive and time-consuming. 
Although, multiphase, multidimensional thermal simulators are appropriate and more 
accurate for the detailed analysis of the process, most of the simulators do not include all of 
these effects. Moreover, the lack of ability to handle sophisticated simulators precludes 
their widespread use in the field. The conventional simulators do not include the 
geomechanical aspects of the process. Absence of flow equations coupled with the stress 
field in the reservoir is the major obstacle for the success of the simulator predictions. The 
most prohibitive factors to use a numerical simulator are cost and time. It may take a few 


months and several hundred thousand dollars to do a comprehensive analysis. 


Before going into a detailed and expensive simulator analysis, a simplistic initial 


assessment of the operation is quite important and valuable. Empirical and correlative 
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methods are mostly specific to a particular field and their predictive ability is also limited. 
Analytical models are capable of rough production estimates for various operating 
conditions. These constitute an excellent tool for process control and operation strategy. A 
reliable analytical model may be useful for the field engineer and can save considerable 
amounts of time and money. Often the lack of available data also makes the use of an 
analytical model imperative. 


A number of analytical models have been proposed in the literature, but most of these 
models are not applicable to CSS in Cold Lake and other heavy oil reservoirs in Alberta, 
where formation parting takes place. McGee, Arthur and Best (1987) and Arthur, Best, Jha 
and Mourits (1991) considered the fracturing phenomenon during the injection period 
which often occurs in oil sands. This was ignored by others who also considered only 
radial flow, whereas the fluid flow from a fractured well is elliptic (Dietrich, 1986). Elliptic 
flow geometry has been investigated in different fields of petroleum engineering: hydraulic 
fracturing, water-flooding, well testing etc. The two extreme views regarding flow during 
production involve the opening and closing of the fracture. Arthur et al. (1991) assumed 
flow into a fracture during production. They introduced anisotropy to handle the opening 
and closure of the fracture but their flow equation was based on an elliptic coordinate 


system which inherently assumes an open fracture. 


In a reservoir, where a considerable amount of primary production exists, the effect of 
temperature change of the incoming fluid from the cold outer zone may not be that 
significant. In the case of heavy oil, it is very important to calculate the heat exchange 
between the hot and cold zone through conduction and convection. Arthur et al. (1991) 
ignored the conduction heat transfer between the two zones. Seba and Perry (1969) also 
considered the calculation of the change of average temperature in the cold zone, but they 


kept the hot zone temperature at steam temperature for the whole project time. 


The most important aspect of the flow calculation is the proper evaluation of oil 
viscosity. So far, there has not been any attempt to calculate an average viscosity based on 
its areal distribution during production. A direct calculation of viscosity from the average 


zone temperature would produce erroneous results. 
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Chapter II 


REVIEW OF THE LITERATURE 


Reduction in the viscosity of the viscous oil around an injecting well is the main reason 
for a higher production response in cyclic steam injection operations. Apart from this 
obvious factor, there are several aspects of the production mechanism which make the 
process complex and diverse. Generation of non-condensable gas, flashing of steam, 
different phase compressibilities, gas evolved from chemical reactions and enhanced 
solution-gas/oil expansion due to increased temperature may contribute to the driving force 
independently or in combination. Solution-gas drive and gravity drainage are two important 
forces in many reservoirs (Prats, 1986; Jones and Cawthon, 1990). Compaction has also 
been cited as the major driving force in Venezuela (de Haan and van Lookeren, 1969). 


The performance of the cyclic steam stimulation process is sensitive to drive and 
drainage mechanisms, reservoir and fluid properties evaluation and operating variables. For 
any predictive model, reservoir temperature calculation and relative permeability 
approximation are particularly important. A small variation in relative permeability may 
significantly alter the production response. Most authors of the several analytical models 
presented in the literature assumed radial flow geometry. Both depletion and non-depletion 
type reservoirs were considered. Many of these works investigated the effect of gravity 
drainage. The role of compaction and solution-gas drive was also explored. 


A few authors used a semi-analytic technique to overcome some of the limitations of a 
fully analytical model. Incorporation of relative permeability and gas phase, non-linear 
variation of fluid heights and transient effects are some of the reasons for taking such an 


approach. 


Steam injection at a thermally efficient rate requires the injection pressure to be higher 
than the formation parting pressure in a variety of Canadian heavy oil and bitumen 
reservoirs. This is particularly true for the earlier cycles. For such cases, analysis of heat 
and fluid flow to and from the fracture are crucial. Although there are a number of papers in 
the literature dealing with heat transfer within the fracture, very few have considered the 


whole reservoir system. 
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It is also important to have a good estimation of the size and shape of the hot or steam 
zone. In many cases this calculation sets the limit of the drainage area. Assumption and 
estimation of the shape of the steam zone are critical for determining the drainage 
mechanism. A number of works in this area are also found in the steam-drive literature. 


Due to the creation of a fracture in some of the heavy oil reservoirs, the flow geometry 
is quite different from that in a radial system. There is a general consensus in the literature 
that the flow geometry around a fractured well is elliptic. Elliptic flow geometry has been 
investigated in different fields of petroleum engineering: hydraulic fracturing, water- 
flooding and well testing. There is a considerable absence of interest regarding its 
application in the cyclic steam injection process. Figure 2.1 gives a schematic overview of 
the work which will be discussed in this Chapter. A quick look at the key features of the 
existing CSS analytical models is presented in Table 2.1. 


2.1 Analytical Models 


The physical models developed in the literature to calculate cyclic steam stimulation 
performance are based on various forces: pressure drawdown, gravity drainage and 
compaction. It is difficult to categorize them in a general fashion. The analytical models are 
reviewed here on the basis of drive and drainage systems. 


2.1.1 Pressure Drawdown 


This approach assumes that there is enough energy in the reservoir to produce some oil 
at the initial reservoir temperature. Apart from reducing oil viscosity in the heated zone, it 
assumes that no drive mechanism is enhanced by steam injection. Discussion of some of 


the major works follows. 


Owens and Suter (1965), in one of the earliest analytical models, recognized the 
inverse proportionality of oil production to the change of viscosity. They assumed that the 
bottomhole temperature followed the harmonic decline trend of the surface temperature. 
They suggested making spot measurements of the bottomhole temperature to calculate the 
parameters of the harmonic equation. Values of temperature and production decline with 
time were presented in a graphical form. Generally, this very simplistic approach tends to 
overestimate the effect of CSS. In a typical case, production enhancement of 20 to 50 times 
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1. Analytical Models 


Pressure Drawdown Gravity Drainage Rock Expansion 


2. Semi-Analytical Models 


3. Steam Zone Growth Calculation 


Frontal Advance Models Override Models 


4. Heat Transfer and Temperature Distribution 


Unfractured Well Fractured Well 


5. Elliptical Flow Regime 


Unaltered Reservoir Composite Reservoir 


Figure 2.1—Review of the Literature 
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would be predicted compared to a model which actually calculates the average reservoir 
temperature. 


Boberg and Lantz (1966) presented the first comprehensive analytical model to calculate 
the production rate of a thermally stimulated well. They found the average temperature of 
the heated zone using correction factors for vertical, radial heat conduction and fluid 
production heat losses. Their radial flow equations are accurate only for conventional 
reservoirs with primary production mainly from the cold zone. Pressure-depleted 
reservoirs, where the production is mainly from the hot zone itself, need a different 
treatment. The heat remaining in the reservoir at the end of each cycle was added to the heat 
injected in the next cycle. They assumed the same initial temperature for the interbedded 
shale, overburden and underburden at the beginning of each cycle. The neglect of low-level 
heating of the oil outside the hot zone and exclusion of any enhanced drive mechanism 
would generally result in an underestimate of their production response. 


Closmann, Ratliff and Truitt (1970) studied the effect of stratification and cross-flow in 
a depletion-type reservoir. The effects of saturation change and viscosity distribution with 
temperature were considered. They neglected gravity and assumed only vertical heat loss 
with a constant heated zone. A Van Everdingen and Hurst (1949) type influx was used to 
calculate the radial flow volume. The way this model is formulated, makes it applicable to a 
first cycle calculation only. 


Ershaghi, Al-Adawiya and Kogawan (1983) presented a graphical method of predicting 
cumulative production. They used the Boberg and Lantz method to generate cumulative 
production data using various parameters. These sets of data were plotted against 
cumulative steam injection on log-log coordinates. It was shown that the slopes of these 
plots were roughly equal to one. All the variables were lumped together into a response 
factor, which was the intercept. There was no rate calculation. Although they matched their 
results with field data, the proportionality of steam injection and oil production cannot be 


applied to many wells. 


Another Boberg and Lantz based model was presented by Gros, Pope and Lake 
(1985). The mathematical formulation of Boberg and Lantz was modified to accommodate 
variable parameters in different layers. They also included temperature dependent relative 
permeability and time dependent average saturation, calculated from a simple material 


balance. No cross-flow between the layers was considered. 
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Because of the need to calculate the size and shape of the heated volume, McGee, 
Arthur and Best (1987) analyzed a vertically fractured reservoir in a steam stimulation 
process. They applied the principle of superposition with Wheeler's (1969) latent heat 
model to find the temperature distribution. A negative value during production and zero 
injection during soak were used. Their main objective was to find out the heated volume. 
For that, they supplied steam injection, volume of oil and water produced, cycle time, 
steam quality and reservoir properties data. The results were used to analyze well 


performance, to calculate well spacing and to determine an overall operating strategy. 


Arthur, Best, Jha and Mourits (1991) extended the work of McGee et al. (1987) and 
presented a comprehensive CSS predictive model for a composite reservoir in elliptical 
flow geometry. They used Boltzman's transformation to solve the flow equations in an 
elliptical coordinate system. The rate equation was an exponential-integral solution for a 
fractured wellbore. They assumed the outer thermal boundary, which was one degree 
higher than the initial reservoir temperature, as their drainage area. They did not consider 
the periodic nature of the 1 coordinate in the elliptic system, thus committing an error in the 
production equations. A more detailed discussion of this work will be presented later. 


2.1.2 Gravity Drainage 


For heavy oil and depletion-type reservoirs, gravity drainage has been identified as a 
critical factor. At the injection period, the steam zone is depleted of oil to some residual oil 
saturation. During the production time, reservoir fluids start flowing into the steam zone 
and resaturate the area. Many authors have employed gravity drainage to calculate this flow 
rate in the heated zone. In heavy oil reservoirs, where it is difficult to move the oil, steam 
heats up the oil around the wellbore; during production, this mobile oil flows into the well 


under the force of gravity. Some of the more important models are discussed below. 


In their study of a depletion-type reservoir, Matthews and Lefkovits (1956) observed 
that at the stripper stage the gas concentration was so low that gravity remained the only 
driving force. By using an average gas-oil free surface height and the pseudo steady-state 
radial flow equation, they derived a gravity drainage equation for a bounded reservoir. 
Later, this equation was widely used by several authors in cyclic steam stimulation 


applications. 
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Doscher (1966), in one of the earlier analyses on CSS process, discussed its 
importance and qualitatively described a model that involved an inner and an outer 
cylindrical segment. He emphasized the significance of the temperature increase calculation 
of the outer segment. It was argued that after the initial production a considerable amount of 
oil would be produced from this section. In his model the hot zone temperature remained at 
the steam temperature during the whole project life. This might create a large error in the 
flow calculation. He proposed using Muskat’s infinite reservoir equation until the pressure 
transient reached the inner zone boundary, and then using the Matthews—Lefkovits (1956) 
bounded reservoir equation. 


Towson and Boberg (1967) used the gravity flow equation derived by Matthews and 
Lefkovits for a bounded reservoir to develop their gravity drainage model. Their thermal 
model was an extension of the previously published Boberg—Lantz model. They compared 
the rate calculated by the gravity drainage equation with the steady-state radial flow 
equation and used the greater of the two values. A step-wise calculation procedure was 
used where the average fluid height of the hot zone was adjusted at each time step. Unlike 
the Boberg—Lantz model, this model gave more realistic results for reservoirs with no 
primary drive. The authors did not consider any effect of temperature on the cold zone 
which always remained at the initial reservoir temperature. 


Seba and Perry (1969) implemented the concept of Doscher (1966) in a concentric 
double cylindrical reservoir with their gravity drainage model. Their model was restricted 
by several assumptions. The radius of the hot zone remained constant assuming no fluid 
displacement during the injection period. Although they calculated the change in 
temperature of the cold zone, they assumed the heated zone temperature remained constant 
at the injected steam temperature. This overestimates production. Their finding of a 
harmonic production decline may not be valid for other reservoirs. 


Jones (1977) reported a depletion-type reservoir model with gravity drainage. He used 
the Boberg and Lantz physical model and the gravity drainage flow equation for his rate 
calculation. A linear water saturation profile with time was assumed. The cumulative effect 
of changing temperature, fluid levels, relative permeability, fluid viscosities and fluid 
saturations was integrated into his calculation. He used scaling factors to calculate water 


and oil flow rates. These factors were determined with significant uncertainty. 


Gontijo and Aziz (1984) improved the Jones model by using a cone shaped steam zone. 
They used Butler, McNab and Lo's (1979) temperature—viscosity relationship for 
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horizontal wells to get the expression for the pseudo-steady state flow potential which was 
a combination of pressure drop and gravity force. A more rigorous heat balance and heat 


loss during soaking was incorporated into this model. They neglected flow from the cold 
zone. 


Sylvester and Chen (1988) presented another model for a pressure depleted reservoir. 
A pseudo-steady state radial flow drainage equation for both water and oil was developed 
in their work. This was derived on the basis of an erroneous equation. A detailed 
discussion of this work is presented in Appendix A. The flow equation contained only the 
fluid height at the wellbore. They considered the steam remaining from the previous cycle 
when calculating the first water saturation value during production. 


Gozde, Chhina and Best (1989) used Gontijo and Aziz's flow equations in their model. 
They also considered cold oil outflow/influx which depends on the pressure difference 
between the interface and the external (outer) boundary. They calculated the pressure 
variation with time at the interface using Ramey's (1970) unsteady-state flow equation for a 
composite radial reservoir. Superposition was imposed in the outer colder region in a 
stepwise calculation. Heat flow from the cap rock to the reservoir during production was 
also considered. A good match with field data was reported. 


Jones (1992) discussed the reasons why cyclic steam stimulation analytical models had 
not been widely accepted by the industry. He argued that the division of the production 
rates of gravity drainage, radial flow models by a scaling factor of 2 to 4 was due to the 
assumption of the entire open wellbore receiving steam. Another factor which had been 
ignored by previous authors was the growth of an air zone at the top of the formation 
which he called the steam thief zone. By taking into account these two observations, he 
constructed a gravity drainage model using existing equations. He reported good history- 
matches with two California heavy-oil projects and also presented some economic 


analyses. 


2.1.3 Rock Expansion 


One model which used an effective fluid compressibility to account for rock expansion, 


which is the most dominant driving force, is reviewed here. 
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De Haan and van Lookeren (1969) simplified the combined effects of all the drive 
mechanisms by representing them by an effective compressibility. Vertical inflow due to 
steam override was treated by finding the effective well radius corresponding to the fracture 
in the steam zone radius at the top of the sand. Then an effective hot zone radius was 
calculated by averaging the steam zone radius and the effective well radius. The authors 
considered the transient effect by assuming semi-steady state flow. They showed that about 
40 percent of the cumulative production in the Tia Juana heavy-oil field, the Bolivar Coast, 
Western Venezuela, was due to compaction drive. 


2.2 Semi-Analytical Models 


Many authors have combined numerical methods with some analytical calculations. 
Most of these semi-analytic approaches were taken to incorporate various complex 
phenomenon of CSS. These works have made considerable contributions to the 


development and understanding different aspects of the steam stimulation process. 


As the temperature of the injected steam increases, it enhances the depletion mechanism 
by making the dissolved gas in the oil less soluble, thus promoting solution gas drive. 
Davidson, Miller and Mueller (1967) used a solution gas drive model to find the volumetric 
rate of oil and water production in the heated zone. They did not consider the fluid 
production heat loss to calculate the average temperature of the reservoir. They found that 
the reservoir response did not depend on the amount of solution gas available within their 
range of gas-oil ratio values. The effect of temperature on oil relative permeability was 
found to be negligible under 400°F, but was significant over 500°F. The limitation imposed 
by the solution gas-oil ratio makes the use of this model very restrictive. 


Martin's (1967) model used only a single sand and is valid primarily for the prediction 
of the first cycle production. He neglected the heat loss through conduction during 
injection, which produced a larger steam zone radius. Radial flow equations for both a 
constant pressure boundary and a no-flow boundary were developed. The effect of 
temperature was incorporated with the inclusion of a temperature dependent oil formation 
volume factor. He solved the production rate and temperature equations numerically using 
an iterative procedure. Controlling the bottom hole pressure value, he analyzed the effect of 
flash steam production on the oil rate and concluded that the oil production almost doubled 


when there was no steam production. 
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Kuo, Shain and Pochas (1970) used a finite differences method to calculate the free 
surface for their gravity drainage model. The heated zone was maintained at a constant 
temperature and the effect of continuous conduction heating of the colder part of the 
reservoir was studied with the help of finite differences. It was concluded that a large 
amount of oil flowed from the area adjacent to the hot zone in the early life despite the 
formation near the outer boundary remaining cold. They presented the results for one cycle 
in terms of cumulative oil production only. The assumption of a constant hot zone 
temperature results in an overestimate of the oil production rate. 


De Swaan (1972) predicted the unsteady three-dimensional temperature distribution 
analytically and used it with the fluid-flow equations for a numerical solution. His model 
could predict the production response of a depletion-type stratified reservoir or of a 
homogeneous gravity drainage reservoir. He noted the importance of radial heat 
conduction. It was concluded that a slight rise of the outer zone oil temperature increased 
the production significantly. 


Williamson, Dake and Chappelear (1976) used an isothermal simulator to do a 
fieldwide investigation of the steam-soak process. They used a semi-analytical approach to 
obtain fluid production rates from the wells by dividing the reservoir into cylindrical 
annuli. Using an average value of the saturation calculation and the flow rate equation, they 
used the Boberg and Lantz method to calculate the average temperature of the hot zone. 
This model was developed to circumvent the difficulties the authors faced when using a 


thermal simulator for a field-wide simulation. 


Flow in the early stages of a CSS cycle is transient in nature. Bidner, Horclo and 
Vampa (1990) argued that Boberg and Lantz's (1966) assumption of steady-state flow 
from the start of production was not correct. They solved numerically the transient 
diffusivity equation in the cold outer zone. A steady-state equation was assumed for the hot 
inner zone. They combined the two equations at the interface through continuity in a 
stepwise manner. Average temperature was calculated according to Boberg and Lantz 
(1966). Three different boundary conditions at the outer drainage radius were considered: 
constant pressure, no-flow and infinite medium. After 95 days of production, no-flow and 
constant pressure conditions produced the minimum and maximum oil respectively, the 
infinite medium assumption gave a result in between these two values. They were not able 
to show any significant difference in cumulative oil production with their approach as 
compared to the Boberg and Lantz method, although early production varied for all three 


boundary conditions. 
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2.3 Steam Zone Growth Calculation 


Production rate and ultimate recovery depend heavily on the size of the steam zone. 
Therefore, it is very important to have a good estimate of the steam or hot zone size. Steam 
zone growth models have been presented in the literature for both cyclic steam stimulation 
and steam-drive processes. They can be classified into one dimensional frontal advance 


models and two dimensional (r-z) steam override models. 


2.3.1 Frontal Advance Models 


In frontal advance models, the steam front is assumed to be vertical in a horizontal 
reservoir. No variation of saturation and temperature in the vertical plane is allowed. For 
horizontal reservoirs, where gravity plays a critical role, these models would overestimate 
the steam zone growth. Review of the principal works is presented below. 


Marx and Langenheim (1959) introduced the first frontal advance steam zone growth 
theory which was widely used in later works on analytical CSS models. A step function 
approximation of the temperature profile for one dimensional radial flow was made. The 
heat balance equation, which considered only vertical heat loss, was solved by Laplace 
transform producing a solution involving the error function. Only latent heat transfer from 


the steam was assumed in their model. 


Ramey (1959) showed that there was no geometric restriction to Marx and 
Langenheim's solution and it could be applied for the variable injection cases with the help 
of superposition, though Ramey’s equation for variable rates was incorrect. In the 
summation term he considered only the superposition of the heat injection terms, ignoring 


the superposition of the time dependent error functions. 


Hearn (1969) modified Marx and Langenheim's solution to include the sensible heat of 
the saturated steam for the heat loss calculation. He concluded that, at a later time for a step 
function temperature profile, enough heat would not be supplied from the latent heat part 
for a limited injection rate. This resulted in a smaller steam zone than previously predicted 


by the Marx—Langenheim solution. 
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Mandl and Volek (1969) examined Marx and Langenheim's limiting assumption of no 
heat flow across the condensing front and introduced the concept of critical time to the 
frontal advance theory. They argued that after the critical time, convective heat transfer 
across the front becomes dominant and cannot be neglected. They presented a steam zone 
growth rate equation after the critical time which depends on reservoir thickness, injection 
rate, temperature and steam quality. 


Closmann (1967) developed the steam zone growth equation for multiple-layer 
reservoirs. The steam zones were assumed to be of equal thickness and heat loss to the 
impermeable layers was also included. In a later work (1968) he studied the growth of 
steam in a preheated reservoir. This was an extension of Marx and Langenheim's work 
which could be employed to calculate the steam zone for the second and later cycles. He 
assumed a linear variation of the initial temperature in the overburden and underburden 
rock, 


Yortsos and Gavalas (1981a, 1981b) developed integral equations for the steam drive 
process and solved them with the help of upper bounds, approximate and asymptotic 
solutions. These methods produced better estimates of steam zones for intermediate and late 
time calculations. They also solved the steam zone growth problem for low quality steam 
and a low injection rate as encountered in the Mandl—Volek model which was valid for a 
convection dominant processes (high injection rate). Calculation of upper bounds on the 
basis of the latent heat ratio of the injected heat is also useful in processes which produce 
receding fronts (such as CSS). The general approach of their model can be used also for 
three-dimensional steam zone growth. 


McGee et al. (1987) calculated the heated zone volume by determining the time at 
which the steam advanced to the formation. This was done by equating the heat injection 
rate and the heat dissipation rate from the fracture face to the formation. They used 
equations derived by Wheeler (1969). Their analysis did not include any flow model and 
needed prior knowledge of production quantities. 


Gajdica, Brigham and Aziz (1990) developed one and two-dimensional semi-analytical 
steam drive models from Yortsos and Gavalas' integral equations. They reported a good 
match with thermal simulator results except when the length to height ratio was less than 


unity. 
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2.3.2 Override Models 


Steam overlay, bypass or gravity override models consider the effect of gravity on the 
density differences of the fluids. Apart from gravity, reservoir layering can also cause an 
override. The relevant papers are presented here. 


Neuman (1975, 1985) proposed a gravity override model where the frontal movement 
is vertically downward. Areal growth of the steam zone was calculated by making a heat 
balance of the heat injected, energy produced and steam condensation rate per unit area. 
This model is applicable only under some specific conditions where the horizontal pressure 
gradient in the steam zone is much less than the vertical pressure gradient caused by the 
difference between the densities of the steam and the liquids. He presented a method to 
calculate the oil displacement beneath the steam zone and the vertical sweep. 


Van Lookeren (1983) studied the effect of gravity override on steam zone growth and 
presented a one-dimensional equation for variation of steam height. An expression for an 
average steam height was derived which could be used to calculate the vertical sweep. He 
indicated that for a very viscous oil (over 1000 mPa.s), viscous tonguing may cause 
considerable bypassing of oil by steam. Thus, his model seems to work well for low and 
medium viscosity oils. There was no attempt to find the volumetric rate of the oil 
production. 


Vogel (1984) studied a descending steam zone model where steam frontal movement is 
vertically downward from the overlaying steam. He proposed a simplification of previous 
models by assuming instantaneous steam overlay for the whole project period. This 
enabled him to use the well known linear heat flow equation in an infinite plane. He 
showed that the error for such a simplification over a ten year period would be only 2.6 


percent for the ultimate heat requirement. 


2.4 Heat Transfer and Temperature Distribution 


Heat transfer calculations are very important for the overall performance prediction of 
the stimulation process. Determination of the temperature distribution in the formation for 
different steam injection schemes at various stages of the process have been studied by 
several authors. These investigations can be classified into two categories—unfractured and 


fractured wells. 
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2.4.1 Unfractured Well 


Because of the complexity of the heat and fluid flow, in most of the studies, many 
approximations and simplifying assumptions were made to find the temperature profile 
around an unfractured well. Several key papers are reviewed here. 


Lauwerier (1955) is credited as being the first to give an approximate analytical solution 
for the temperature distribution in the reservoir and the surrounding formation. He assumed 
a linear one-dimensional reservoir with constant physical and fluid properties. The 
temperature distribution within the reservoir was independent of vertical position and heat 
flow within the reservoir was by convection only. He only considered vertical conduction 
heat loss. 


Boberg and Lantz (1966) in their model considered temperature distribution during the 
production period. Heat loss was calculated for vertical and radial conduction and by 
production of hot fluid. The temperature was represented by a step function where the 
unheated zone remained at the initial reservoir temperature. They assumed the reservoir 
sand and adjacent rock heat capacities to be the same. An equation for average reservoir 
temperature was presented which used approximate correction factors for all three modes of 
heat loss. This simple, easy and fairly accurate method attracted several authors who used it 
in their steam stimulation analytical models. 


Davidson et al. (1967) developed an equation for the average temperature in the heated 
region during the soaking period. For the production period they only considered vertical 
and radial heat loss. These quantities were approximated by some equivalent expressions 


and a dimensionless temperature distribution in differential form was presented. 


De Haan and van Lookeren (1969) also assumed conduction heat losses in the vertical 
direction only and used Carslaw and Jaeger's (1959) simple analytical solution to calculate 


the average temperature of the steam zone. 


Closmann et al. (1970) derived a simplified expression for the temperature distribution 
in the reservoir after steam injection by assuming only vertical conduction heat loss. They 
assumed the thermal properties of oil bearing sand and the adjacent cap and base rock to be 


equal. 
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2.4.2 Fractured Well 


Thermal recovery in many heavy oil reservoirs requires high steam injectivity. In very 
viscous oil-bearing formations, it can be achieved only by parting the formation. 
Especially for CSS, heat transfer around a fractured well and calculation of the temperature 
profile during the injection and soaking period is critical for performance prediction. A 
review of the more important papers is presented here. 


Dysart and Whitsitt (1967) first attempted to calculate the transient heat transfer from a 
vertical fracture by making an energy and mass balance of an elementary volume of the 
fracture. Heat flow to the overburden was neglected and only conduction heat transfer to 
the pay sand was assumed. In a subsequent paper (Whitsitt and Dysart, 1970) the authors 
improved the model by allowing fluid leak-off from the fracture face. They assumed a 
linear variation of the leak-off from zero at the wellbore to the maximum value at the 
fracture tip. This model was a representative of a fracture which used additives to minimize 
leak-off. Heat loss from the top and bottom of the element was neglected again. Despite the 
over-simplification of the model, it did serve the purpose of identifying some aspects of 
hydraulic fracturing. 


Wheeler (1969) presented a comprehensive analytical investigation of heat transfer from 
a fracture. He studied three different models—hot water injection, steam injection and a 
latent heat model. The hot water model accounted for convection in the fracture, conduction 
and convection in the pay zone and conduction in the overburden. A steam injection model 
also considered condensation in the vertical fracture. He predicted the presence of a non- 
isothermal water zone in front of an isothermal steam zone in the fracture. This model is 
valid for a thick impermeable formation since heat losses to the adjacent formation are 
neglected. His latent heat model did not allow for a hot water front in the fracture but heat 
loss to the overburden was considered. Equating the injected heat rate to the energy 
transferred to the sand by conduction and convection, he neglected any energy 
accumulation in the fracture. This assumption creates a condensate advancing front from 
the fracture face instead of a steam advancing front. He argued that the error due to this 
limitation would probably be no worse than assuming a step function temperature 
distribution in the fracture. A constant injection rate was also assumed. Wheeler reported a 
close match of his steam injection model results with a more sophisticated numerical model. 
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Sinclair (1971) used Wheeler's latent heat model which assumed a constant leak-off 


velocity from the fracture face. He studied the effects of heat transfer in deep well 
fracturing. 


Ben-Naceur and Stephenson (1985) reviewed different heat transfer mechanisms 
during fracture propagation for their numerical model. They considered temperature- 
dependent fluid properties in the model. They cited a cubic approximation for temperature 
penetration in the formation which matched fairly well with the results from the analytical 
and numerical methods. 


Biot, Massé and Medlin (1987) determined that, during hydraulic fracture growth, the 
fracture fluid temperature linearly rises to the reservoir temperature until reaching the 
middle of the fracture and then remained constant beyond that point. They used two 
different methods to obtain the solution: the method of characteristics and the method of 
successive approximation. 


Meyer (1987) presented an analytical solution which combined transient heat 
conduction and convection from the fracture with a convection boundary condition at the 
fracture face. His model also included a time-dependent fracture face temperature and 
energy storage. He coupled the energy and two different fracture propagation models (i.e. 
Geertsma-de Klerk and Perkins-Kern-Nordgren) to find a closed-form dimensionless 
integral solution. This medel seems to be the most realistic in handling the actual heat and 
fluid flow mechanisms in the fracture. 


Boone and Bharatha (1993) studied the temperature distribution along a propagating 
leak-off dominated fracture. They developed a temperature profile along a horizontal 
fracture assuming only conduction heat loss. Later they extended this to find two limiting 
bounded solutions which included convectional heat loss. Their solutions are based on the 
assumptions of Carter’s (1957) leak-off dominated fracture and Marx—Langenheim's steam 
injection theory. They concluded that, in Cold Lake oil sands, effective permeability is the 
primary factor for the advancement of heat. The initial heat transfer is dominated by 


conduction and the effect of convection is felt by the later cycles only. 
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2.5 Elliptical Flow Regime 


Flow around a fractured well and in anisotropic reservoirs has generated interest in 
elliptic flow equations. Studies of the elliptical flow regime have been made for 


homogeneous and composite reservoirs. A brief description of the key papers is presented 
below. 


2.5.1 Homogeneous Reservoir 


The earliest flow studies in the elliptic coordinate system were for homogeneous 
unaltered reservoirs. Some of these papers are reviewed here. 


Coats, Tek and Katz (1959) investigated the unsteady state problem of water 
encroachment into elliptic shaped reservoirs. A numerical solution of the diffusivity 
equation in the elliptic coordinate system was presented. They showed that the water influx 
calculation was about 7% in error when an elliptic boundary was approximated by an equal 
area circle for short dimensionless time. 


Prats (1961) and Prats, Hazebroek and Strickler (1962) considered the elliptic flow of 
incompressible and compressible fluid respectively in a reservoir with an infinite- 
conductivity vertical fracture. For an incompressible fluid, the choice of an elliptic 
boundary for a large area (req>2) has a negligible effect on the pressure transient response. 


In the later work it was shown that the effect on compressible fluid could also be ignored. 


A closed-form solution for flow into a circular wellbore in an elliptic drainage system 
was proposed by van der Ploeg, Kirkham and Boast (1971). Their work was related to 
water flow in a confined elliptical aquifer. A free surface was assumed at the water head. 
Steady-state solutions were developed for various well locations using gravity flow. 


Results and flow nets were presented for several cases. 


In 1974, Gringarten, Ramey and Raghavan studied the unsteady-state pressure 
distribution created by an infinite-conductivity vertical fracture. This study was done 
because of a lack of predictive models to analyze the performance of short-time pressure 
transients. They applied Green's function to obtain both early-time and long-time solutions. 
In comparison with a radial system, they concluded that the effective well radius for an 
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infinite conductivity vertical fracture in an infinite reservoir during the pseudoradial (long- 
time) period is exactly one-fourth of the fracture length. 


Kucuk and Brigham (1979) reported a study where they tried to handle elliptical 
problems in a unified and consistent manner. They gave analytical closed-form solutions 
for constant pressure and constant rate wellbore conditions. A long time approximation was 
also presented. In their following work (1981) they applied the elliptic flow solution to a 
water influx problem. For a constant pressure interface they constructed tables and curves 
of dimensionless influx and rate for different ag/be ratios. Using different effective radii for 
a closed-form solution, they concluded that for intermediate and long times, an effective 
radius based on elliptic flow gave the most accurate results. 


Hale and Evers (1981) used a conformal transform to solve the steady-state elliptic 
problem for application in vertically fractured gas wells. By defining a radius of 
investigation which changes with time, they modified the steady-state solution to 
approximate the unsteady-state solution for both the constant pressure and the constant rate 
case. 


Perkins and Gonzales (1985) investigated the effect of thermoelastic stress on injection 
well fracturing. They discussed the change of stress due to the temperature difference of 
injected fluid and the reservoir fluid. They used Marx and Langenheim's method to 
calculate the volume of the elliptic shaped flooded zone and derived an expression to find 
the major and minor axes of the ellipse. 


2.5.2 Composite Reservoir 


Water and steamflood processes (or any other drive method) create two distinctly 
different regions in the reservoir with different fluid properties. Investigations of these 


composite reservoirs are briefly discussed below. 


Obut and Ertekin (1987) developed a transient flow equation for a composite system in 
elliptic coordinates. Assuming equal pressure and flux at the interface they gave solutions 
for both constant pressure and constant rate at the wellbore using separation of variables. 
The well is an infinite-conductivity vertical fracture. Their solution is a good approximation 


of the stated problem. 
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The same problem was addressed by Stanislav, Easwaran and Kokal (1987). They 
gave a solution that is strictly applicable for equal diffusivity in both regions, but they 
showed that even for a diffusivity ratio as high as four, they got fairly good results. They 
attributed this to the more robust stability properties of Mathieu's function with respect to 
changes in parameters. In the case of steam stimulation, where the value of diffusivity ratio 


may be much higher than four, their equation fails to produce any meaningful result. 


One study on dual porosity reservoirs was reported which followed the same approach 
as that used for a composite reservoir. Okoye, Matese, Hayatdavoudi and Ghalambor 
(1988) gave a solution for an elliptical flow system in dual porosity reservoirs. They 
concluded that elliptic geometry could be more appropriate for natural and vertically 
fractured reservoirs where anisotropy creates an elliptical boundary. 
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Chapter ITI 


STATEMENT OF THE PROBLEM 


With the steady depletion of Canada's conventional oil resources, ever-increasing 
attention is being given to the heavy oil reservoirs. The largest commercial heavy oil 
operation is in Cold Lake, Alberta. Cyclic Steam Stimulation (CSS) is the principal mode 


of thermal operation currently employed in Cold Lake and other heavy oil reservoirs. 


Analytical models are excellent tools to make a quick assessment of new projects as 
well as current operations before a large undertaking. Although the need for such tools is 
well recognized, few models adequately address the problems involved in Cold Lake CSS 
operations. The objective of this research is to explore heat and fluid flow around an open 
fracture, as well as after it closes, and develop an analytical model which incorporates the 
actual flow situation in the field. As oil viscosity is the single most important parameter in 
flow calculations, special attention is given to the proper evaluation of average temperature 
and oil viscosity. 


3.1 Statement of the Problem 


The object of this study is to develop an analytical model which includes formation 
fracturing. For realistic modelling, the following tasks have to be accomplished: 


1. Examination and evaluation of flow patterns around a fracture. 
2. Appropriate heat balance and flow calculations. 

3. Proper evaluation of oil viscosity in reservoir. 

4, Computer implementation of the solution. 


After a program is written to implement the ideas, the model should be able to 
successfully predict Cold Lake and other heavy oil CSS performances. It should give 
proper flow patterns and accurate temperature (and oil viscosity) values. It is also expected 


to appropriately respond to changes in a wide variety of parameters. 
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Chapter IV 


MODEL DEVELOPMENT 


One of the most significant applications of the Cyclic Steam Stimulation (CSS) process 
is in very heavy oil reservoirs. There are a number of analytical models reported in the 
literature which deal with lighter oils. All of these models consider radial flow into a well. 


In the case of heavy oil reservoirs like Cold Lake, Alberta, the flow conditions are quite 
different from conventional reservoirs. Steam injection rates under normal conditions are 
too low and injection pressures must be increased until parting occurs to achieve a high 
injection rate. Fluid flow, heat transfer and mass balance have to account for the presence 
of this fracture at the wellbore. There are very few works which deal with this process in 
the CSS literature. McGee et al. (1987) probed the importance of a fractured wellbore in a 
fieldwide CSS operation. Arthur et al. (1991) presented the only analytical predictive 
model. 


The mechanisms of fluid and heat flow in CSS are complex and diverse. Several 
authors (Vittoratos, 1991; Denbina, Boberg and Rotter, 1991; Beattie, Boberg and McNab, 
1991) have discussed different aspects of this complex process. Although chemical 
reaction, relative permeability hysteresis, compaction and fluid expansion with temperature 
are present in CSS to various extents, gravity drainage and depletion type reservoirs have 
generated the most interest for their significance in heavy, as well as light oil reservoirs. 


Although the importance of gravity drainage into a fractured well has been recognized, 
so far no one has attempted to investigate it. Elliptic flow occurs in the area surrounding a 
vertical fracture in an anisotropic formation or in an aquifer with an elliptic inner boundary 
(Kucuk and Brigham, 1979). In the case of pressure drawdown, Arthur ef al. (1991) 
assumed the existence of a vertical fracture in an isotropic formation during the injection 
period, yet they solved the flow equation assuming an anisotropic reservoir for elliptic 
flow. Although anisotropy was introduced for the fracturing phenomenon, directional 
permeability was not included in the flow equation. To close the fracture during 
production, Arthur ef a/. assumed an isotropic reservoir, essentially treating the reservoir 
as an elliptical system with an open fracture at the wellbore. In reality, the reservoir remains 
isotropic for most of the time. The effect of any directional permeability would be 
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significant during production and limited to an area close to the fracture face. The case of a 
closed fracture should also be examined as a part of an elliptic drainage system. 


4.1 Flow Equations 


Three different flow equations are considered in the present study. The vertical fracture 
is assumed to be of negligible width and is bounded by over- and underburden. A 
conceptual diagram of the fluid and heat flow is shown in Figure 4.1. The reservoir is 
divided into two concentric elliptical cylinders. The inner region is called the hot zone and 
the outer region is called the warm zone. All three geometric configurations for the flow 
equations are presented in Figure 4.2. 


4.1.1 Flow into a Wellbore at the Centre of an Elliptical Drainage Boundary 


It may be assumed that the fracture formed during injection is completely closed during 
production without leaving any residual permeability in the hot zone. Under this 
assumption, the flow geometry in the reservoir is an elliptical shaped confined hot zone 
with a well at the centre (Figure 4.2a). Van der Ploeg, Kirkham and Boast (1971) 
developed a closed-form solution for steady saturated flow into a fully penetrating well in 
an elliptical flow geometry. They considered several well locations for their study. The 
essence of their approach was to derive orothonormal functions for the problem. The 
coefficients Anm were evaluated for various cases. If one were concerned about flow at the 
wellbore only, a few values of ANm were needed for quick convergence. These were used 
to determine Ano, which was the only coefficient required in the flow equation. A well 


located at the centre with reservoir thickness of h has a flow rate of, 
_ —20KhA®PAy, 
of = 
Ps 
where A® is the difference in hydraulic head, K is the hydraulic conductivity of the 
reservoir, a is the major axis of the ellipse and rw is the wellbore radius. A detailed 


description of the development and solution procedure of this equation is given in 
Appendix B. 


(4.1) 
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Cold Zone 


Plan view 


Heat & Fluid Flow 


Heat Flow 


Front view 


Figure 4.1— Conceptual diagram of fluid flow and heat transfer 
for the proposed model. 
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a. Circular well at the centre of an elliptical drainage system 


Yo 


b. Approximation of the elliptical drainage system with a well at the centre 


Fracture Face 


c. Linear gravity flow approximation with a fractured wellbore 


Figure 4.2—Plan views of different flow geometries. 
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4.1.2 Flow From a Finite Line Source into a Well 


If an approximation of the above equation were to be made, one could imagine a long 
rectangular block producing oil from a well at the centre. As the system is symmetric, only 
one-half of it needs to be considered. Muskat's (1937) equation for flow from a finite line 


source can be used (Figure 4.2b). For a well, which is on the perpendicular bisector of the 
line source, the flow equation is, 


2nkh(p, — P,) 


ee ee eg 4.2 
r c? 


w 


where yo is the distance between the well and the line source and 2c is the length of the line 
source. If (yo/c) << 1, the above equation can be written as, 


ie 2nkh( p, - P,, ) 


win 2% 
ie 


w 


(4.3) 


4.1.3 Gravity Drainage 


From numerical simulation (Dietrich, 1986; Pethrick, Sennhauser and Harding, 1988) 
and observation wells (Duerksen, Cruikshank and Wasserman, 1984; Reis, 1990), it has 
been observed that the heated elliptic zone around the fracture is quite eccentric. An 
eccentricity (€) value of zero represents a circle whereas unit eccentricity is a straight line. 
The perpendicular distance between the boundary of this zone and the fracture face at the 
steam injection point varies from 2 to 8 meters over the life of a well, whereas the fracture 
length could be several hundred meters. So, after injecting steam, essentially a very long 
rectangular block around the fracture with hot mobile oil would be created (Nzekwu, 
Hallam and Williams, 1990). Figure 4.2c shows this flow geometry. Due to the high 
permeability in the fractured zone, linear gravity flow of this hot oil into the fracture would 
take place. Gringarten et al. (1974) reported linear flow from a vertical fracture at early 


times. This occurs close to the fracture face. 


Muskat (1937) described linear gravity flow in conjunction with his discussion of flow 


through dams. He presented a simple equation for the flux per unit length: 
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Forchheimer (1930) proposed this equation for a linear system based on the 
assumptions of Dupuit. This equation is known as the Dupuit—Forchheimer equation. 
Muskat showed that the assumptions made by Dupuit were in marked contradiction to those 
given by calculations with rigorous hydrodynamic theory. On the other hand, he showed 
that the same equation could be derived using an approximate potential theory. He found 
that the results were remarkably close approximations of the values of the flux determined 
using rigorous analytical or empirical methods. 


The above equation is not an exact solution of the problem but is very much within the 
error bounds of the calculation of a rigorous solution. 


4.1.4 Outer Warm Zone 


Steady-state elliptical flow was assumed to be taking place in the outer warm zone at 
each time step. The equation (Muskat, 1937) for this situation can be written as 


_ 2nkh(p, ~ Ps) 


(4.5) 
un( 24>] 
a +0, 


4.2 Pressure 


For an approximate solution to the flow from an elliptical drainage boundary into a well 
at the centre, an exponential decline curve was used in the present study to control the 
wellbore fluid height. The index of this equation was used as a history-matching parameter. 
In the proposed linear gravity models, the wellbore fluid heights were fixed for each cycle. 
These heights were allowed to decline at the rate of 1/(total cycle production time) to 
achieve stability and an improved history-match. Average fluid heights at any time step n 
for the hot and warm zones were estimated from the following equations which were 


derived for the present study 


Nhot — hijos an (Qnor ~ (Tae) / PA,,,.> (4.6) 
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Fears = hwarm — Oss / CAs, (4.7) 
where Q is the total fluid flow at each time step. 


The interface fluid height was estimated from an area weighted average of fluid heights 
in the two zones. This height determined the pressure at the interface. 
fe norA ar Ten A 


hot warm ; (4.8) 


h, 
Ayo + Al i 


Fluid height at the drainage boundary was taken to be constant and hence the pressure 
was assumed to be the initial reservoir pressure. Jones (1977) and Chen (1987) reported 
very little change in fluid height at the interface for a smaller radial system. Considering 
low flow rates and the huge size of the outer warm zone, this assumption would not 
produce any appreciable error. The change in the interface fluid height was monitored to 
track any gross error in this assumption. 


4.3 Fracture Length 


If no fluid flows from the fracture into the formation, it can be assumed that the fracture 
always extends beyond the calculated heat front in the fracture. For a constant leak-off 
velocity, a simple material balance will determine the maximum extent of the fracture in the 
formation ( Wheeler, 1969). Assuming that only condensate flows into the formation and 


that there is no accumulation in the fracture, the fracture length is given by 


ip 
=—#-, (4.9) 
: AVP. 
This equation, in effect, assumes that the fracture tip always propagates ahead of the 
heat front in the fracture. Although there would be a time when this condition will not be 
satisfied, that time is usually large. Wheeler reported achieving accurate results for injecting 


steam at 159 m3/day for 100 days. 


The other option which has been included in the model is Carter's (1957) fracture 


propagation equation, which was derived as a function of time assuming constant width, 
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where C is the fluid loss coefficient, & is effective fracture width and Sp is the spurt loss 
quotient. For large injection times, the above equation reduces to the familiar form of 


LINE, 


Pare (4.11) 


For shorter injection time or small fluid leak-off, Carter's equation overestimates the 
fracture length. It gives a fair estimate of the fracture length for most of the cyclic steam 
stimulation projects. There is provision in the model to specify fracture length for each 
cycle. 


4.4 Initial Temperature Calculation 


Steam is injected into the formation to decrease the oil viscosity. Although there are 
several analytical models which deal with the temperature profile inside the fracture, very 
few include fracture and overburden together in the calculation. Wheeler (1969) presented 
the only comprehensive analytical calculations for heat transfer from a vertical fracture. He 
presented three models based on different states of fluids and modes of heat transfer. The 


following assumptions were made for the "latent heat" model. 


1. Fracture height is equal to the sand height. 
2. Two equal and symmetric fracture wings of negligible width are formed. 
3. The pay sand is infinite horizontally and confined by impermeable overburden 
and underburden strata vertically. 
4. Mode of heat transfer: 
a. Fracture - convection and condensation. 
b. Sand - conduction and convection. 
c. Cap and base rock - conduction. 
5. No energy accumulation in the fracture. 
6. Temperature in the fracture is defined by a unit step function. 
7. Leak-off is constant for the injection period. 
8. Energy from the fracture into the formation is transferred by condensate. 
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Given the injection rate and total injection time, the following equation gives the 
temperature profile for any point in the reservoir including the fracture, pay sand and over- 
and underburden. 


\Blv*(r-7,)|+— 


(Roel 
oo| oe Qv./T | is (4.12) 


This equation is used to determine how temperature is distributed in the hot and warm 
zones after the injection period. The dimensionless time for the steam front to reach any 
point in the fracture was estimated by 


We 2Hs\ Tc) for [Pe flo dv. (4.13) 


The temperature in the fracture maintained a step function profile. Initially only 
condensate flows into the formation. Once the energy supplied to the fracture exceeds the 
energy transferred into the formation, a steam zone will advance. The time for this process 


was Calculated from the following equation, 


2X;Afp 
ae aa 


(4.14) 


we ierfo(a/T,, ) = 


Using these three equations presented by Wheeler (1969), an isotherm of 100°C was 
used as the hot zone boundary in the present study. The warm zone boundary was selected 
to be the isotherm at 0.01 degrees above the initial reservoir temperature. The choice of 
100°C is based on the flowing mobility of the oil at Wolf Lake as reported by Pethrick et 
al. (1988). In general, this temperature may be defined on the basis of oil mobility for a 
particular reservoir. A sensitivity study of the choice of boundary temperature will be 


presented later. 


The average temperature of the two zones was estimated using a linear polynomial 


approximation (trapezoidal). 
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4.5 Heat Loss Equation 


Carslaw and Jaeger (1959) presented a transient solution for a heat dissipating semi- 
infinite slab inside an infinite medium. The slab had a unit initial temperature. They solved 
the one dimensional heat conduction equation in the vertical direction to obtain the solution 
which involved the error function. Boberg and Lantz (1966) used this solution as part of 
their heat loss calculation to determine the average steam zone temperature during 
production. Marx and Langenheim (1959) used a similar approach to find the steam zone 
growth during injection. 


Souza, Pedrosa Jr. and Marchesin (1987) used an analytical-numerical method to 
determine the heat loss to the overburden and underburden. This method was developed for 
an easier treatment of the heat loss equation in the overburden. They employed an 
approximation to solve the heat conduction equation in a finite medium with variable 
temperature at the reservoir-overburden boundary. This eliminated the need for an 
additional grid system for the overlying strata and storage of temperature history of each 
grid block. 


Vinsome and Westerveld (1980), with an aim to save computer storage and grid block 
calculation, adopted a simpler strategy. They argued that due to the uncertainty in the cap 
and base rock thermal properties, high accuracy in the solution of the heat loss calculation 
was not required. They proposed a much simpler semi-analytical solution of the linear one- 
dimensional conduction equation based on the results found from the more complex 
variational method (Weinstein, 1972, 1974; Chase and O'Dell, 1973). The assumptions for 


their solution were, 


1. Temperature in the cap rock varies smoothly. 
Horizontal heat conduction in the cap rock can usually be neglected. 
3. The tail of the temperature distribution contains little energy and is relatively 


unimportant. 


Based on the above assumptions, they gave a fitting function for the temperature profile 
which adequately met all the boundary conditions. Inserting this function into the one- 
dimensional heat conduction equation at the interface using a finite difference discretization 


of the time derivative, the following equation was found 


= N 
Bate Wg, 4-424) (4.15) 
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where @, is the thermal diffusivity, To is the interface temperature, d is the diffusivity 
length and p and q are parameters. The parameter ToN is the interface temperature at the 
beginning of the time count. This ensured that the heat loss rate decreased from cycle to 
cycle in a CSS implementation of this method. This also eliminated the need of 


superposition. Heat loss rate and energy stored in the cap and base rock were given by, 


T. 
qu = k=) (4.16) 
E, = “ d(f + pd +2qd’) (4.17) 


h 


where k, is the thermal conductivity. This solution was used in the present model in order 
to calculate the heat loss to the over- and underburden as well as the heat conduction 


between the two zones. Appendix C presents the equations which were used in evaluating 
the average temperatures. 


4.6 Average Soak and Production Temperature 


Vogel (1984) presented a simple heat-loss equation for the over- and underburden 
assuming that the steam covered the entire heated area instantly. Several authors (Sylvester 
and Chen, 1988; Gozde et al., 1989) have used this equation to calculate the heat loss 
during the soaking period. 


t 

Hogs: = 2k, 77; (Ts -— Tp )| 2 4.18) 

los h At S y) = ( 

A cylindrically shaped heated zone was assumed in developing the equation. This 
quantity of heat was deducted from the total heat injected to determine the heat available at 
the beginning of a production period. Although, the rate of heat loss is very high during the 
early stages, very few authors tried to find the change of average temperature during the 


soaking period. 


The most popular method for calculating the average reservoir temperature during 
production is the Boberg and Lantz technique. Performing a heat balance for a cylindrical 


heated zone, they derived the following equation, 
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Tiyge = Ty + (Ts ~ Ty)[9-6:(1 - 5) — 5] (4.19) 


where, 9, and @; are the unit solutions for conduction in the radial and vertical directions, 
respectively, and 6 gives the fraction of the heat loss due to fluid production. The values of 
6, and @, were presented graphically against dimensionless time. The solution given for 
@, was valid for dimensionless time greater than 0.1. Bentsen and Donahue (1969) 
proposed an alternative solution in terms of a series function which was applicable for a 
dimensionless time of less than 0.1. They also presented a unit solution for one- 
dimensional vertical heat loss through conduction. The parameter 6 was found by 
integrating the heat loss due to fluid production at each time step. Gontijo and Aziz (1984) 
produced the unit solutions with the help of regression analysis. This was done for 
programming simplicity. 


Arthur et al. (1991) treated the problem differently. They included a second zone, just 
outside the hot zone, in their overall heat balance. The outermost boundary temperature was 
taken to be one degree above the initial reservoir temperature. By so doing, they eliminated 
the need for calculating the radial heat loss. The average temperature of the warm (outer) 
zone is important for fluid flow calculation in the later part of a cycle. Although they 
considered heat exchange between the zones through fluid flow, they ignored conductive 
heat transfer. The inclusion of conduction is especially important during the soak period 


when no fluid movement takes place. 


A new set of equations was developed (Appendix D) in the present study which 
included heat conduction between the two zones. All conductive heat transfer terms were 
calculated using the approximation suggested by Vinsome and Westerveld (1980). A more 
detailed discussion of their work is presented in Section 4.5. The heat loss rates were 
evaluated with time-dependent surface temperatures. The need for superposition was 


avoided by using cumulative project time. 


Following are the overall heat balances to calculate the average zone temperatures 


during the production and soak periods. 


Hot zone heat balance: 
Heat content= —_ production heat loss + over- and underburden heat loss 


+ warm zone heat loss — heat gain from warm zone flow. 
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Warm zone heat balance: 


Heat content = heat loss to hot zone flow + over- and underburden heat loss 
— heat gain from hot zone. 


Heat exchanges due to fluid flows during the soak period were considered to be zero. 


4.7 Saturation 


Earlier CSS analytical models did not involve any fluid saturations, because single 
phase flow was assumed. Most of these models evaluated only the oil production rate. 
Water production was calculated from an arbitrary water-oil ratio. No consideration was 
given to saturation changes or relative permeability. 


Jones (1977) introduced an analytical method to evaluate saturation changes based on a 
water material balance in the hot zone. He correlated these saturation values with relative 
permeabilities measured in the laboratory. Following Jones, other authors (Gontijo and 
Aziz, 1984; Sylvester and Chen, 1988; Gozde et al., 1989) used his model or a modified 
version of it. 


Composite elliptical flow models (Stanislav et al., 1987; Obut and Ertekin, 1987), 
developed mainly for weli-testing, used single-phase flow. Hence, the need for saturation 
calculation did not arise. Arthur et al. (1991) employed two different zones to do an overall 
material balance in both oil and water to calculate the saturation at each time step. Based on 
their mobile fluid calculation, two different schemes were developed and tested in the 


present study. 


In the first scheme, the traditional approach taken by most authors was followed and an 
overall material balance of only water in both the hot and warm zones was considered. 
Although there is evidence of the presence of a gas phase in many reservoirs, the amount of 
gas in heavy oil reservoirs may be small. It was also assumed that most of the injected 
steam condensed during the soak period. Both zones were represented by average zone 
values. Initial mobile fluid saturations were calculated from the following general equation. 


Initial mobile fluid = initial volume + injected volume + cum. inflow 
—cum. production — cum. outflow. (4.20) 
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If the fluid saturations in each zone did not add up to unity, the values were normalized. 


Subsequent saturation values were calculated from a simple material balance on the water, 


emir boul Gy not Iu 1 Fina)! OV po) (4.21) 
ete es (oer Ott G5) (OV a) (4.22) 
ae ey (4.23) 
Ss =1—S..4 (4.24) 


where Swm1 and Swm2 are the initial mobile water saturations in the hot and warm zones, 
respectively. 


For the second scheme, a material balance on both oil and water was found. It was 
assumed that the total production time was much larger than the individual production time 
step. Saturation equations were developed from the continuity equations assuming the 
change in hot zone volume to be equal to the change in warm zone volume. Arthur et al. 
(1991) normalized the saturation values at each time step, whereas more weight was given 
to the water balance in the present study. The oil saturation was calculated from the 
saturation equation at each time step instead of normalizing the values. Arthur et al.'s 
mobile saturation equations and the modified calculation methods for this study are 
presented in Appendix E. For both schemes, the hot zone water saturation was restricted to 
1 - Sor. Mobile water, consequently, was restricted to a maximum value of 1 - Sor - Swe. 


4.8 Average Viscosity Calculation 


Viscosity is the single most important parameter in the calculation of thermally 
stimulated oil production calculation. Boberg and Lantz (1966) developed their flow 
equation solely on the effect of change in oil viscosity in the hot steam zone. They 
incorporated this effect by increasing the effective wellbore radius. An earlier model by 
Owens and Suter (1965) used a viscosity-temperature graph to estimate the oil production 
rate. Closmann et al. (1970) evaluated the time average viscosity values to use in a van 
Everdingen and Hurst (1949) type oil influx calculation, which involved superposition. 
They did not address the areal distribution of viscosity. Since these early models, all 
authors have used the average temperature of the hot zone to calculate oil viscosity. The 


reason for following this procedure is the lack of any simple analytical equation to calculate 
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a temperature distribution during the production period. Viscosity values found from a 
viscosity-temperature relationship corresponding to the average zone temperatures cause an 
overestimation of oil production. These average values do not give the proper weight to the 


effect of lower temperatures at the outer edge of the averaging area and hence, give lower 
viscosity values. 


_ Ay + Aly + Ast, + Alt, 


Have = 
A, (4.25) 


(A, = A, + A, + A, + A,) 


An indirect, simple but effective approach was taken in the present study to overcome 
this limitation. Using Wheeler's (1969) equation, it is possible to calculate all the isotherms 
at the end of an injection period. As the strips between these isotherms are linear beds in 
series with the same cross-sectional area, the average viscosity corresponding to the 
average temperature of each strip was calculated. True area-weighted averages of the oil 
and water viscosities for each zone were found corresponding to different overall average 
temperatures. This was accomplished by varying the steam injection time, rate and steam 
quality. To develop a correlation between the average temperatures and the average 
viscosities, the total area was kept constant for all calculations. Third degree polynomial 
equations were derived from the generated data for both the water and the oil. The 
temperature range covered by these equations was adequate for all the variation in 


temperature during the production time. 


4.9 Flow Reduction Factor 


Using the full values of length and height of the fracture in the flow equation produces 
unrealistic flow rates. In theory, one would get such high rates, but in the actual reservoir 
there are several factors which would reduce the flow drastically. Three such factors are 


discussed below. 
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4.9.1 Vertical Steam Sweep Efficiency 


Due to gravity override, the injected steam does not sweep the whole pay thickness. 
This is true irrespective of the creation of a fracture during injection. Myhill and Stegemeier 
(1978) examined the oil-steam ratios predicted by different frontal displacement models in 
steamflood calculations. They compared the results with actual field data and laboratory 
experiments. The measured values were generally less than the predicted values. They cited 
steam override in the real situation, which was not addressed by the models they studied, to 
be one of the factors for such lower values. Van Lookeren (1983) presented several 
calculations and field data to corroborate this. Even without a thief zone, the sweep 
efficiency may be as low as 20%. Very seldom does it exceed 70%. Jones (1992) 
discussed the reasons why analytical cyclic steam stimulation models have not been 
accepted widely. He showed with the help of a temperature survey in an observation well 
near an in-progress steam cycle well that only a fraction of the open pay zone took steam. 
He argued that the common assumption of the entire open wellbore receiving steam 
injection was one of the major reasons why the predictive flow equations overestimated the 
flow by factors of two to four. Reis (1990) reported that only 40% of the completed 
interval accepted steam in a California heavy oil reservoir. 


4.9.2 Fracture Parameters 


The complexities of steam injection as well as reservoir heterogeneity make the 
prediction of fracture parameters unique for each oil sand reservoir. Interpretation and 
analysis of field data play a significant role in characterizing a fracture. Single or multiple 
well operation, temperature, viscosity, grain size, thief zone, depth of burial, clay content 
and chemical reaction are among the many other factors contributing to the uniqueness of 
each oil sand reservoir. The most significant fracture parameters are the height and the half 
length. Although these parameters may be altered significantly by localized geological 
effects, general trends in Alberta heavy oil are discussed below. 


4.9.2.1 Fracture Height 


Arthur et al. (1991) presented the only CSS analytical model which included a vertical 
fracture at the wellbore. They used the total pay thickness (100%) for the fracture height. 
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Although this ideal situation may exist in some reservoirs, several authors reported partial 
vertical climb of fracture propagation in Alberta heavy oil. 


Leshchyshyn, Faroug Ali, Settari and Chan (1994) reported a fracture height of less 
than 50% of continuous net pay height. This was attributed to shale heterogeneity. In 
another work, Leshchyshyn and Kennedy (1994) inferred from a post core clay analysis 
that the vertical steam fracture climbed at an angle of 45 degrees. This was also supported 
by an independent analysis of Formation Micro-Scanner (FMS) logs. Lowe and Donovan 
(1989) observed a vertical fracture profile around the wellbore and steeply rising as it 
moved out from the wellbore in Athabasca oil sand. Fialka, McClanahan, Robb and 
Longstaffe (1993) analyzed a post steam core from a Cold Lake area oil sand reservoir. 
After four cycles, they found the vertical sweep efficiency to be only 60%. 


4.9.2.2 Fracture Length 


In a fracture growth analysis, Leshchyshyn, Seyer and Langdon (1991) found that only 
60% of the calculated fracture length was active with hot fluid (of temperature above 90°C). 
Many authors assume a complete closure of the fracture during production. Reis (1990) 
reported formation of vertical fractures in a deeper, less permeable and more consolidated 
heavy oil reservoir in California, which remained open during production. Chan and 
Sarioglu (1992), in their discussion on the complexities of the CSS process, noted that 
numerical simulation studies for different reservoirs suggested the presence of closed but 
not completely healed fracture. Moreover, fracture length calculated from an analytical 
relationship is only approximate. In many cases, field observations showed two to five 
times overestimation of the fracture length calculated by such an analytical formula 
(Leshchyshyn, 1995). 


Flow reduction factor values between 0.20 to 0.50 were used in the present study. This 
incorporated the effects of all the factors discussed above. These values, along with 


wellbore fluid heights, were used as fitting parameters to history match the results. 
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4.10 Relative Permeability 


Only two-phase flow was considered in this study: oil and water. This assumption is 
reasonable when most of the injected steam condenses before the start of the production 
period and the gas liberated in a heavy oil reservoir is considered negligible. 


It is not possible to find a correlation of water and oil relative permeabilities with water 
saturation which would be applicable to a wide range of reservoirs. Even with the help of 
experimental data for the same reservoir, it is difficult to match field production data. 
Gontijo and Aziz (1984) used Farouq Ali's (1982) correlation for a Midway-Sunset field 
which was originally developed for a Kern River steamflood (Gomma, 1980). Chen and 
Sylvester (1990) showed that the values found in the Gontijo and Aziz (1984) study were 
considerably different from the reported field data and these correlations caused significant 
errors in the oil rate calculation. 


In addition, large discrepancies may occur due to hysteresis. Bang (1984) was only 
able to match field production by adjusting relative permeability in numerical simulation 
studies. He found that the relative permeability to water at residual oil for imbibition 
(injection) was about seven times greater than that for drainage (production). Dietrich 
(1981) also reported that relative permeability curves generated through numerical 
simulation to reproduce observed water-oil ratios in several CSS operations worldwide, 
were much lower than those of routine laboratory imbibition water-oil relative permeability 


tests. 


The effect of temperature on the water relative permeability end point is well 
established. Faroug Ali (1982) concluded from numerical simulation of (Cold Lake, 
Alberta) viscous oil that the temperature effect on water-oil relative permeability must be 
taken into account for meaningful history-matching. Dietrich (1981) and Nakornthap and 
Evans (1986) presented numerical simulation results showing an increase in the irreducible 
water saturation with temperature. This shift in end point causes a decrease in the water 


relative permeability values. 


For simulation of the Wolf Lake project, Pethrick et al. (1988) adjusted relative 
permeability values, found from experiments, within the bound of the errors. For the 
present study, the tabulated values of water-oil relative permeability presented in their work 
were used. These values were adjusted for change in water relative permeability end point 


due to the effect of temperature. 
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4] 
Chapter V 


SOLUTION AND COMPUTER IMPLEMENTATION 


A program was written in FORTRAN to implement the ideas described in the previous 
chapter. The main purpose of this program is to calculate flow rates at each time step during 
production. Average temperatures of the hot and warm zones at the end of the injection and 
soak periods are also evaluated. Special attention was paid to avoid iterative procedures. 
No “trial and error” approach was taken at any level of the calculation. Empirical formulas 
and correlations were used to minimize data requirements. Whenever needed, cubic spline 


interpolation and a table look up routine was used to find intermediate values of tabulated 


data. 


5.1 Description of the Model 


ib 


Saturated steam is injected into the wellbore at high pressures, causing formation 
parting. Fracture half-length is calculated either from a given fluid leak-off 
coefficient (Carter, 1957) or material balance (Wheeler, 1969). 


The reservoir is divided into two concentric elliptical cylinders around a vertical 
fracture at the wellbore. The inner cylinder, which extends from the fracture face 
to the 100°C isotherm, is called the hot zone. The outer cylinder, which extends 
from the heated boundary to the drainage boundary, is called the warm zone. 
The drainage boundary is the isotherm just above the initial reservoir temperature 
(16.0 +0.01°C). 


At the end of injection, average temperatures and the distances of the isotherms 


separating the zones are determined. 


During soak and production, conduction heat loss to the over- and underburden is 
calculated. There is no radial heat loss from the outermost boundary. Average 
zone temperatures are calculated from energy balance equations considering heat 
loss due to fluid production, heat loss to the confining strata, and convection and 


conduction heat exchange between the zones (Appendix D). 
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5. Three different flow equations are used for the three cases to calculate the flow in 
the hot zone. Interface fluid height and the corresponding pressure are calculated 
from a material balance equation. It is assumed that steady-state flow takes place in 
both zones during a time step. Pressure at the wellbore is maintained at a constant 
level during production for the gravity flow models. Although no flow takes 
place across the outer boundary, the fluid height is assumed constant (initial 
reservoir pressure). The elliptical steady-state flow equation is used to determine 
the rate in the warm zone. Despite steady-state flow during each time step, the 
overall flow is transient owing to the pressure variation at the hot-warm interface. 
Flow rates at the beginning of a time step are assumed constant for that step. 


6. Heat remaining in the reservoir at the end of each cycle is added to the injected heat 
in the next cycle, allowing the reservoir temperature to be kept at the initial value at 
the beginning of each cycle. 


5.2 Description of the Subroutines 


The subroutines implementing the model are described here in a sequential manner 


along with the principal features of the main program. 


5.2.1 Initialization 
Subroutine DATAIN 


This is the first subroutine called by the main program. Data are required to be supplied 
in SI units. There are five sets of input data. The first set needs injection, production, and 
soak times for each cycle along with steam pressure, saturation temperature and a fluid loss 
coefficient (optional). The second set is comprised of porosity, formation thickness, 
injection rate, steam quality, absolute reservoir permeability and initial reservoir 
temperature. Compressibilities, saturations and initial reservoir pressure are specified in the 
third set. The fourth set provides thermal conductivities and heat capacities of the reservoir 
and the confining strata. The final set includes all the relative permeability data. Data are 


provided in free format. 


42 


ced gies y th eet thee ae ae hag, 


+ | ks Sen aOR Se tlre 
nha 
bra ntey gon OSGi a7 tedtiew HEP ys cage 
aed were Get bie inaani eaneteal a 

dig © 1.5 (29o ore vital tual soll al aT 
4 Dent bs pte Hehe ve thd ery Beng Leste 

cient haere ep A a ae 

cab ts ar ed nl eer y vier catty Vit 08 Geerney 


vite Tarts Wie tral ee eA Osh B= Fae 


wt Sti Sali We ea 5 Gt they ory 7 Re 


. (ele SANt ney’! hand soit 4H a) ori Oj _ gees ae 
- — 
vue wate fo a 


+"5 


series? 2 to 


c ; gir? si Lote ily Pp iorwapm 
Pi jae OY } . Ciera 


ab OR ols a pines wiunetiehttie 

Si iw te oT es aes 
‘A. css cee “eae REE (OF 

sdegttty ice hy mans ep Al Lgeba os: dgee 

Jiw oe. (UL eR Seine gaia. ag 


yc ce Albis way a ae 


Wh tee aT hy a 


Subroutine SPLINE 


This subroutine calculates all the spline coefficients for the tabulated relative 
permeability data. 


Subroutine OVISC 


Oil and water viscosities and densities are calculated with the help of this program. 
Viscosity values at initial reservoir temperature are calculated here before the start of the 
injection period. 


5.2.2 Injection Period 


The main program initializes all the project and cumulative values at this stage. 
Maximum and minimum production limits of 100 m3/day and 1m3/day, respectively, are 
also set at this time. The maximum number of cycles is supplied and all the cycle 
calculation values are initialized. 


Subroutine PVTIN 


This routine calculates the viscosity, density and heat capacity of the steam, water and 
oil. The latent heat of vaporization is also evaluated. At this point, all the values are 


calculated at steam temperature. 
Subroutine TEMPGAUSS 


Hot and warm zone boundaries are calculated by this program, along with the fracture 
length. The fluid leak-off velocity can be found in three different ways: supplying of leak- 
off coefficients as direct input, calculating the leak-off velocity through the subroutine 
LEAKOFF or using Wheeler’s (1969) material balance equation with a known fracture 
length. Wheeler’s equations are integrated using 16-point Gaussian quadrature. Average 
temperatures at the end of injection for each zone are also calculated here. In the case of an 
unknown fracture length, it can be found either from Carter’s or Wheeler’s equation. For 
repeat cycle calculations, the heat remaining from the previous cycle is added in the form of 


a dimensionless quantity to the heat input in the present cycle. 
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Subroutine PRESSURE 


Pressures at the interface and the outer boundary are calculated for the approximate 
solution to the flow into a wellbore at the centre of an elliptical drainage system. Pressure at 
the wellbore is controlled by an exponential function of production time. 


5.2.3 Soak Period 
Subroutine SOAKTEMP 


This program evaluates the average zone temperatures at the end of the soak period. 
Parameters for the fitting function are first evaluated here. It also calculates the heat loss to 
the over- and underburden during this period. A time step size of 1 day was used for the 
short soak period. Moreover, at high temperature the heat loss rate varies considerably for a 
small change in temperature. A 1 day time step enabled retention of better accuracy. 


At this point , the main program calculates all the initial mobile saturations and PVT 
data at average temperatures. Saturations for the first production time step are evaluated and 
limits are checked. The oil in place in the hot zone is also calculated. Then the production 
stage calculation values are initialized. 


5.2.4 Production Period 


Relative permeabilities and mobilities for the first production time step are evaluated 
here. In addition to the table look-up (subroutine SEARCH), relative permeability 
calculations from correlations are provided for. Wellbore fluid height and all flow rates are 
calculated next. Maximum and minimum flow rates are checked. Cycle cumulative values 


are updated. 
Subroutine SATN 


This program calculates the saturation values based on the flow rates evaluated at the 
present time. Options for an oil- and water-based material balance calculation as well as 
only water based calculation are provided. The new saturation values are used for the next 


time step. 
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Production results for the time step are written here. Then fluid heights at the hot and 
warm zones as well as at the interface are updated. 


Subroutine AVPRODTEMP 


This routine evaluates the average zone temperatures at each production time step. The 
conduction, heat flow from the hot zone to the warm zone is multiplied by a scale-down 
factor. This is done in accordance with the assumptions made for the flow reduction factor. 
Fifty percent of the formation thickness was considered active in the heat exchange. Heat 
lost to the over- and underburden and through fluid production is also calculated by this 
program. 


The main program is looped with a time increment of two days until the end of 
production. Cycle results are then written and project cumulative values are updated. The 
heat remaining in the reservoir is calculated at this time before the start of the next cycle 
calculation. The procedure is repeated until the end of the cycles is reached. A flow diagram 
of the computer program is shown in Figure 5.1. 
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Figure 5.1—Computational flow diagram of the model. 
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Chapter VI 


RESULTS AND DISCUSSION 


This study examines Cyclic Steam Stimulation (CSS) performances under fractured as 
well as unfractured wellbore conditions. Flow into a completely closed fracture which 
involves Elliptic Flow into a circular Well (EFW) is investigated. A closed-form as well as 
an approximate solution to this problem is examined. Based on the results of these studies, 
a partially open fracture situation is investigated with the help of the new model discussed 
in Chapter IV. To verify the effectiveness, as well as the limitations of the new model, an 
extensive sensitivity study is conducted. The effects of changes in downhole steam quality, 
formation thickness, steam injection rate, time step size, steam soak period, steam 
temperature, hot zone boundary temperature and initial reservoir temperature are discussed. 
Two different methods of calculating the saturations are examined. The effect of 
temperature change on the water relative permeability end point is also studied. 


The new model investigates the effects of conventional method of evaluating the oil and 
water viscosities at each time step during production. Results from a new and much 
improved method, as discussed in Chapter IV, are presented. The effect on the average 
temperature due to conduction heat loss from the hot inner zone to the warm outer zone is 


investigated. 


Sample results are discussed and interpreted for each of the above mentioned aspects of 
the CSS process in this chapter. Figure 6.1 shows the investigation structure of this study. 
The objectives are to explore the existing concepts and to provide a better understanding of 
the CSS process. The relative importance of different phenomena on the overall CSS 


process is investigated. 


6.1 Simulation Data 


The most frustrating part of this study was the lack of available field data. All the past 
and present operators in Cold Lake, Wolf Lake and Marguerite Lake heavy oil pilot and 
commercial projects were contacted. Data were also requested from AOSTRA and ERCB. 


All these efforts came to naught. Either the data were confidential or they were on sale fora 
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Simulation Runs 


Closed fracture Fractured wellbore 
Wellbore 
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Elliptical Flow Elliptical Flow Linear Linear 
into a Well into a Well Gravity with Gravity with 
(Approx.) (Closed-form) Total Balance Water Balance 
(EGE) (LGWB) 


Base Run (LGTB 


Adjusted Rel Perm.) 
New Average Viscosity Heat Conduction Betn. Zones Sesitivity Runs 


Hot Zone Boundary Temperature 
Soak Period 
Increasing Steam Volume 
Slug Size 
Time Step 
Steam Temperature 
Steam Quality 
Initial Reservoir Temperature 
Formation Thickness 


Figure 6.1—Overall scheme of the simulation runs conducted. 
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price beyond the means of the present undertaking. This study was conducted using the 
only published data (Pethrick er al., 1988) on the Wolf Lake project, jointly owned by 
Petro-Canada and BP Resources Canada Limited. The simulation input data are shown in 
Table 6.1 and the oil-water relative permeabilities are shown in Table 6.2. The correlations 


used to calculate the fluid properties are presented in Appendix F. 


The fracture equations do not account for geological stresses or change in stress. The 
criterion for the fracture model is same as specified by Pethrick ef al.(1988) that the 
pressure at the leading edge of the fracture is 8500+500 kPa at the end of the injection 
period. No injection restriction was imposed. This is the approximate fracture propagation 
pressure. 


At BP’s former Marguerite Lake lease, bitumen is produced at an approximate depth of 
450m (Buckles, 1979) primarily from the lower Cretaceous Clearwater formation. Wells 
are predominantly in 220 mx 100 m five-spot configuration. Based on the net present value 
of the project, Pethrick et al. optimized the injection and production strategies of individual 
wells in Wolf Lake projects. Prior to this optimization, several years of CSS production 
data from BP’s Phase A Pilot in the Wolf Lake project were matched to validate the thermal 
simulator. The agreement between the simulation results and field behaviour allowed a high 
degree of confidence in the optimization strategy. The final strategy in a 220 m spacing is 
shown in Table 6.3. This data has been used as the base run for all the experiments in the 
present work. 


For the initial investigation and the sensitivity study, the first five cycle results are used. 
Several production performance parameters are used to analyze the results. Project 
efficiency is compared on the basis of the oil-steam ratio (OSR). The operational criteria are 
based on a calendar day oil rate (CDOR). The remaining efficiency parameter used is the 
water-oil ratio (WOR). 


6.2 Elliptical Flow into a Well (EFW) 


Two different approaches were taken to study this problem. The approximate method 
was tried to verify the merit of further analysis by the closed-form solution. The closed- 


form solution required some more rigorous calculation. 
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Table 6.1—Simulation input data. 


Porosity, fraction 0.305 
Permeability, m2 1.0 
Initial reservoir pressure, kPa 2750.0 
Initial oil saturation, fraction 0.66 
Initial reservoir temperature, °C 16.0 
Irreducible water saturation, fraction 0.34 
Residual oil saturation, fraction 0.35 
Rock compressibility, kPa-! 0.45E-6 
Water compressibility, kPa7! 0.45E-6 
Oil compressibility, kPa™! 0.34E-6 
Oil thermal expansion coefficient, °C-! 0.848E-3 
Reservoir specific heat, kJ m3K-! 2350.0 
Reservoir thermal conductivity, kJ m!d4!K-! 149.6 
Cap/base rock specific heat, kJ m3K-! 2350.0 
Cap/base rock thermal conductivity, 

kJ m-!d-1K-1 149.6 

Oil viscosity at 16°C, mPa.s 8 1000.0 

Oil viscosity at 100°C, mPa.s 115.0 

Oil viscosity at 200°C, mPa.s 6.8 

Net pay, m 23.8 


Fracture half length, m 240.0 
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Table 6.2—Water-oil relative permeability (Pethrick et al., 1988). 


0.3200 0.0000 1.0000 
0.4300 0.0025 0.6440 
0.5450 0.0100 0.3910 
0.5929 0.0200 0.2420 
9.635957 0.033027 021380 
0.6760 0.0530 0.0610 
0.7145 0.0730 0.0300 
0.7530 0.1130 0.0120 
0.7700 0.1200 0.0000 


Table 6.3—Operational data. 


Steam Slug Injection Soak/down Production Cycle Time 
Cycle Size (m3) Time (days) Time (days) Time (days) (days) 


1 6250 ZS 15 89 129 
7000 28 15 123 166 
7500 30 15 135 180 
8000 32 15 150 i 


8500 34 15 174 BL 


zZ 
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6.2.1 EFW Approximate Solution 


Initial investigation of this problem was done with the help of an approximate solution. 
An exponential decline curve, Pwr = Pj exp(-t/d), was used to calculate the bottomhole 
flowing pressure. The value of Py was fixed at 1500 kPa below the steam injection 
pressure for each cycle and d was varied from 260 in the Ist cycle to 90 in the 5th cycle for 
best matching. The interface pressure of this model was approximated using an infinite 
integral solution and a material balance equation. This was based on Arthur et al.’s (1991) 
solution procedure. 


Results of the approximate as well as the closed-form solution are presented in Figure 
6.2 through Figure 6.8. The analysis for the closed-form solution will be presented in the 
next section. Figure 6.2 shows the cycle oil production along with BP's simulation results. 
A fairly good match was obtained by adjusting the wellbore fluid pressure. Predicted 
production in the 3rd cycle is about 100 m3 less than the simulation value. An over- 
production of about the same amount in the 4th cycle kept the cumulative oil production 
very close to the simulation result. This can be seen in Figure 6.4, where both the BP 
simulation and the approximate solution predict the same total amount of oil at the end of 
five cycles. In the history-matching effort, more emphasis was given to match the oil 
production and let the water production be according to the relative permeability 
distribution. Figure 6.3 shows the cycle water production of the EFW schemes and the BP 
simulation run. It can be seen that the predicted water production by the approximate 
solution is consistently below the simulation values for all the cycles. Although the general 
trend of increasing water production is seen, no attempt was made to adjust the relative 
permeability curves. Any adjustment for hysteresis or temperature effect would have 
produced even less water. Cumulative water production at the end of five cycles is only 
9,042 m3 compared to 23,724 m3 projected by the BP simulation (Fig. 6.5). A new set of 
relative permeability values could be generated to force a match. The objective of the study 
was to remain as close to the original simulation input data as possible and explore the 


actual results of various schemes. No forcing of the solution was done for this reason. 


Because of the incorrect estimation of water production, one cannot put much 
confidence in the performance parameters predicted by this model. Cycle OSRs (Fig. 6.6) 
and CDORs (Fig. 6.7) matched the values predicted by the BP the simulation fairly closely. 
The effect of low water production can be seen in Figure 6.8, where the WOR varied from 
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Figure 6.2—Comparison of cycle oil production 
for closed fracture models. 
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Figure 6.3—Comparison of cycle water production 
for closed fracture models. 
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Figure 6.4—Comparison of cumulative oil production 
for closed fracture models. 
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Figure 6.5—Comparison of cumulative water production 
for closed fracture models. 
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Figure 6.6—Comparison of cycle oil-steam ratio for closed 


fracture models. 
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Figure 6.7—Comparison of calendar day oil production rate 
for closed fracture models. 
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Figure 6.8—Comparison of cycle water-oil ratio 
for closed fracture models. 
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a low of 0.5 in the Ist cycle to about 2.0 in the Sth cycle. These values are somewhat lower 
than the reported average field results (Faroug Ali, 1994) from the Cold Lake area. 


6.2.2 EFW Closed-form Solution 


The approximate method gives an inconclusive answer regarding flow geometry. 
Although the oil production was matched satisfactorily, water production was far below the 
simulation results. This discrepancy did not support a correct interpretation of the flow 
geometry assumption, as all input data were unmodified. Before getting into any 
mechanistic analysis, it was decided to conduct further investigation with the help of a 
closed-form solution of the problem (van der Ploeg et al., 1971). 


The closed-form solution requires the computation of the coefficient ANg. This 
involved integration of several functions which are described in Appendix B. A Gaussian 
quadrature method was employed to do these integrations. The calculated geometric 
parameters of the inner hot ellipse were used as the input data (a=240 m, b=2-3.5 m, 
Tw=0.2 m). The average value of ANo was evaluated to be 2.05. After observing low flow 
rates in all the cycles, the wellbore fluid height was fixed at 0.0 m to achieve the maximum 
flow potential. 


Figure 6.2 compares the cycle oil production of both the EFW models with the BP 
simulation. Although the shape of the closed-form solution curve is similar to the 
approximate solution, the flow volumes are the lowest for the closed-form solution. While 
the water production for the 1st cycle is as high as the simulation result, it drops drastically 
in the 2nd cycle (Fig. 6.3) and falls below the approximate solution after that. The 
cumulative oil production is only half of that predicted by the BP simulation (Fig. 6.4). It 
should be noted that the approximate solution uses pressure drawdown as the driving 
force. To achieve a higher production rate to match the BP simulation results, three 
different options were available to calculate the average pressure of the hot zone: average 
pressure calculated from an infinite integral solution (Arthur et al., 1991), saturated 
pressure corresponding to the average temperature and average pressure corresponding to 
the average fluid height of the hot zone. The maximum of these values was used. 
Moreover, higher flow reduction factors were used to match the results. Tamim and Farouq 
Ali (1995) reported less oil production also for the approximate model using the same flow 
reduction factors for both models. This produced even lower rates for the closed-form 
solution. The total water production at the end of the Sth cycle is substantially higher than 
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the approximate solution but about half of the simulation results (Fig. 6.5). Higher 
prediction by the closed-form solution was due to the large amount of water production in 
the Ist cycle. 


Due to the low oil production, both cycle OSRs (Fig. 6.6) and cycle CDORs (Fig. 6.7) 
were much lower than the simulation values. The 3rd cycle produced the lowest CDOR of 
4.4 m3/day as well as the lowest OSR of 0.107. Interestingly, the WOR remained constant 
at 2.0 for all cycles (Fig. 6.8). The reason for such behaviour is not understood. This is 
close to the average WOR for the first 5 cycles in other Cold Lake area operations (Farouq 
Ali, 1994). 


Lower production rates for the EFW models indicate that a circular well at the centre of 
an elliptical drainage boundary does not represent the true flow geometry in the field. Any 
fracture at the wellbore which is created during injection does not completely close during 
production. Even if it is closed, it does not heal and considerable high residual permeability 
remains in the fracture plane. Although discrepancies in oil and water production volumes 
are evident, curves for all the parameters followed the same general trend for both cases. 
The cycle water production and consequently the WOR curve followed different directions. 
This could be attributed to the large difference in 1st cycle water production between the 
closed-form and the approximate solutions. Average temperature and saturation calculations 
for a subsequent cycle depend on heat remaining in the reservoir and the last saturation 
values from the previous cycle. The relative permeability curves and the nature of the flow 


equation determine the Ist cycle response. 


6.3 New Model 


From the investigation of flow in an elliptical drainage system with a circular well at the 
centre, it is found that the flow rates are much lower than the BP simulation results. 
Although the approximate solution predicted satisfactory oil production, this was achieved 
using much higher flow reduction factors than the probable field values. Both the 
approximate and the closed-form solution predicted much lower water production than the 


BP simulation results. 


A new flow geometry is proposed for the new model. It is assumed that the vertical 
fracture created during the injection period does not close completely during the production 
period or it leaves behind a vertical plane with a high residual permeability. Linear gravity 
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flow will take place from a long narrow rectangular block of hot mobile oil into this highly 
permeable vertical plane. For the nature of the flow equation, this model will also be called 
the Linear Gravity flow model. The flow reduction factor was increased from 0.25 in the 


Ist cycle to 0.50 in the 8th cycle. The wellbore fluid height was fixed at 19.5 m for all 
cycles. 


6.3.1 Saturation Schemes 


Two different saturation calculation schemes have been tested. In the first scheme, in 
accordance with the technique used in all radial flow analytical models, oil and water 
saturations during production are calculated based on an overall material balance of water in 
both the hot and the warm zones. This model is called the Linear Gravity flow based on 
Water Balance (LGWB). For the second scheme, the material balance was done on both the 
oil and the water. This will be called the Linear Gravity flow based on Total Balance 
(LGTB). 


Cycle oil productions for the two saturation models are shown in Figure 6.9. Model 
LGWB consistently produced slightly more oil at each cycle. Model LGTB produced the 
same amount of oil for each of the first three cycles, after which the oil production fell off. 
Overall, both models were able to predict the oil production fairly closely to the results 
given by the BP simulation. The maximum difference in production was only 7%. This can 
be seen more readily in Figure 6.10 showing the cumulative oil production. The production 
forecasted by the LGWB model was only 200 m3 more than the BP simulation after 5 
cycles. The accuracy for water production was not as good. After the 2nd cycle, the water 
production kept declining for the LGWB model. On the other hand, the LGTB model 
produced a considerably more water than the BP simulation (Fig. 6.11). It overestimated 
water production by 125% in the first cycle and about 60% in the 5th cycle. It should be 
mentioned that no modification of the relative permeability curves was made for either 
model. Emphasis was placed on matching the oil production and history-matching 
parameters were adjusted independently for each model. Because of the higher production 
in the 1st cycle, the cumulative water production for model LGWB was within a reasonable 
range of the BP simulation value whereas, model LGTB predicted 60% more total water 
production. Figure 6.12 displays the cumulative water production for all the models. 


The production performance parameters found from these models cannot be accepted 
for making operational decisions. Close matching of oil production generated good 
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Figure 6.9—Comparison of cycle oil production 
for fractured wellbore models. 
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Figure 6.10—Comparison of cumulative oil production 
for fractured wellbore models. 
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Figure 6.11—Comparison of cycle water production 
for fractured wellbore models. 


4 10° 
—Q— BP Simulation 


3.5104 | — -ModelLGTB 


— © - Model LGWB 
3 104 


2.5 104 
2 10° 
1.5 104 


1 10* 


3 
Cycle No. 


Figure 6.12—Comparison of cumulative water production 
for fractured wellbore models. 
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agreement of the related parameters — OSRs (Fig. 6.13) and CDORs (Fig. 6.14) but, the 
WOR for LGTB was two times higher than the BP simulation due to much higher water 
production. Interestingly Figure 6.15 shows that LGWB has a constant WOR of 2.0 for all 
the cycles. The EFW closed-form solution also produced the same WORs. No explanation 
could be found for such a behaviour. 


It is evident from the results that just considering the water balance of the system in 
calculating saturation does not produce the expected amount of water. The constant value of 
2.0 is contrary to the increasing WOR with increasing cycle numbers as reported in almost 
all field observations as well as in the BP simulation. On the other hand, model LGTB 
produced from one and a half times to two times more water for each cycle. Both the 
cumulative water production and the WOR curves look like they have been shifted upward 
from the BP simulation values. Neither of these saturation calculation schemes meets all the 
boundary conditions as rigorously as in a numerical simulation. Each model gives 
approximate values. Moreover, they are also somewhat sensitive to the time step size. More 
accurate values are found from smaller steps, although larger steps do give acceptable 
results (detailed analysis is presented in the sensitivity study). 


6.3.2 Effect of Temperature on Relative Permeability 


Frizzell (1990) did an extensive analysis of 15 years of thermal laboratory data to 
determine the effect of temperature on relative permeability and the saturation end-points in 
heavy oil reservoirs. He cited several studies and concluded that the irreducible water 
saturation increased with increasing temperature. Observing the sub-parallel nature of the 
WOR curve of model LGTB and the BP simulation in Figure 6.15, intuitively one would 
conclude that incorporating the effect of temperature on the water end-point saturation 
would improve the results. Several adjustments were made to include this effect in the 
relative permeability curve and Figure 6.16 shows the final form of the adjusted relative 
permeability values. These were not changed or readjusted for the entire project time. A 
cubic spline interpolation scheme was used to determine the intermediate values. 


Results from this experiment are presented in Figure 6.17 through Figure 6.23. Oil 
production from each cycle is shown in Figure 6.17. It can be observed that a very good 
match was achieved with the BP simulation. Excellent agreement can also be seen in the 
cumulative oil production plot (Fig. 6.18). Figure 6.19 shows the cycle water production 
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Figure 6.13—Comparison of cycle oil-steam ratio 
for fractured wellbore models. 
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Figure 6.14—Comparison of calendar day oil production rate 
for fractured wellbore models. 
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Figure 6.15—Comparison of cycle water-oil ratio 
for fractured wellbore models. 
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Figure 6.16— Adjusted relative permeability. 
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Figure 6.17— Oil production at each cycle for model LGTB, 
with adjusted relative permeability due to 
temperature effect. 
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Figure 6.18—Cumulative oil production for model LGTB, 


with adjusted relative permeability due to 
temperature effect. 
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from model LGTB with and without the effect of temperature on the water end point 
saturation. Improvement of the water production forecast can be seen readily. A fairly good 
match is achieved up to the Sth cycle. The adjusted LGTB still produced a much larger 
volume of water for the last three cycles. Such high water production could be due to 
relative permeability hysteresis. However, it is believed that the hysteresis effect is more 
pronounced in the early cycles, whereas the water production is higher in the last three 
cycles. Two possible explanations for later cycle high water production are found in the 
literature. Vittoratos (1990, 1991) attributed this to the formation of a well defined single 
phase emulsion flow in which the water content rises with cycle number. Other authors 
also reported a similar observation (Gallant, Stark and Taylor, 1993). Karyampudi (1995) 
reported the change of reservoir rock wettability near the wellbore from water-wet to oil- 
wet, because of exposure to high temperature steam. This causes reduced oil mobility and a 
higher water cut. Adjustment to the relative permeability curve reduced the cumulative 
water production by 30% up to the Sth cycle (Fig. 6.20) . 


Excellent agreement was achieved for both the OSR (Fig. 6.21) and the CDOR (Fig. 
6.22). For the first four cycles, WORs are approximately equal to the actual cycle numbers. 
Figure 6.23 shows the WOR. Excessively high WORs are seen for cycles 7 and 8. They 
can be disregarded for any operational or economic analysis. Figure 6.24 shows the oil 
recovered as a percentage of the original oil in place. At the end of 8 cycles, which is 
towards the end of the operational life of a well in Cold Lake, 15% of the oil has been 
recovered. This is close to the average ultimate recovery of 20% in Cold Lake projects 
(Boone, Gallant and Kry, 1993). Very little oil was produced from the warm zone in all 


runs. 


The overall performance of the adjusted LGTB model was very good. It successfully 
matched most of the BP simulation results. The first five cycles of this model have been 


used as a base run for all the sensitivity runs. 


6.4 Effect of Viscosity 


The most significant parameter for flow calculation in thermal oil recovery processes is 
the oil viscosity. So far, no one has attempted to calculate an average viscosity based on the 
actual areal distribution of temperature during production. Results from the two methods of 
calculating viscosities are presented in this section. In the direct method, the average zone 


temperature is found from a heat balance equation at each time step and the corresponding 
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Figure 6.19—Water if beet at each cycle for model LGTB, 
with adjusted relative permeability due to 
temperature effect. 
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Figure 6.20—Cumulative water production for model LGTB, 
with adjusted relative permeability due to 
temperature effect. 
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Figure 6.21—Oil-steam ratio for model LGTB, with adjusted 
relative permeability due to temperature effect. 
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Figure 6.22—Calendar day oil production rate for model 
LGTB, with adjusted relative permeability 
due to temperature effect. 
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Figure 6.23— Water-oil ratio for model LGTB, with adjusted 
relative permeability due to temperature effect. 
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Figure 6.24— Predicted oil recovery as a percentage of 
Original oil in place. 
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viscosity is calculated from a viscosity-temperature relationship (i.e. =aeb/T). The average 
method uses a viscosity-temperature relationship based on an area weighted temperature 
distribution for the same average zone temperature. A detailed description of the method is 
presented in Section 4.8 of Chapter IV. 


Average temperatures of both the hot and the cold zones are calculated using Eqn. 4.12. 
Isotherms are found by increasing the value of the y coordinate (€ in dimensionless form), 
which is the perpendicular distance from the fracture face. Increase of € is continued until 
the 100°C isotherm is reached, and then the average hot zone temperature is calculated 
using the trapezoidal rule. The step size of this incremental distance was varied to check the 
sensitivity of the average temperature (average viscosity). For the 1st cycle, dimensionless 
step sizes of 0.02 (0.476 m), 0.01 (0.238 m), 0.008 (0.1904 m) and 0.005 (0.119m) 
yielded average hot zone temperatures of 190.62°C, 197.49°C, 198.41°C and 199.0°C 
respectively. When the dimensionless step size of 0.005 needed 17 steps to reach the 100°C 
isotherm, step size 0.02 needed only 5 steps to reach the same isotherm. Decreasing the 
step size beyond 0.005 (0.119 m) did not change the average temperature much. The 100°C 
isotherm varies from 2 m in the Ist cycle to about 5 m in the 8th cycle . Step sizes may 
gradually be increased with the cycle number without loosing accuracy. Slightly higher step 
size (0.015) may change the cyclic production responses but it would not cause any 
considerable error in ultimate recovery (less than 5%). Step size of 0.005 has been used for 
the 1st cycle in the present study. 


The oil viscosity calculated by the direct method is much lower than that from the 
average method. At 134°C the direct method estimates an oil viscosity 6 times lower than 
the average method. At higher temperatures, this difference is reduced considerably. Figure 
6.25 shows the variation of oil viscosities with temperature calculated by both methods. It 
can be seen that even at 178°C the direct method underestimates the oil viscosity value by 3 
times. The reason for this discrepancy is the lack of any provision in the direct method to 
accommodate the influence of a rapidly decreasing viscosity values towards the hot zone 
boundary. At higher temperatures, this difference is smaller than at lower temperatures. 
This is due to the highly non-linear nature of the viscosity-temperature relationship. The 
average method gives a proper weight to the higher viscosity values near the hot zone 
boundary in calculating the overall zone viscosity. Clearly, this method produces more 
accurate and representative results. Because of the almost linear nature of the water 
viscosity-temperature relationship, differences in the water viscosity remain essentially 
constant across the temperature range. Figure 6.26 shows only 30% underestimation in the 
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Figure 6.25—Comparison of area weighted average 
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water viscosity calculation by the direct method. This is due to the small variation in water 
viscosity with change in temperature. The consequence of lower oil viscosity on 
computation by the direct method can be seen in Figure 6.27 where the cumulative oil 
production has been plotted against time. After 673 days the direct method produces 3,689 
m3 of extra oil which is 45% more than the average method. Most of this oil is produced in 
the early cycles when the maximum oil is available in the hot zone (Fig. 6.28). 


A higher production in early cycles portrays a deceptively high efficiency. Figure 6.29 
shows the cycle OSRs for both the methods. For the first two cycles, the direct method 
shows around two times more oil production for the same amount of steam injected in the 
average method. The CDORs are also much higher for the first two cycles. Figure 6.30 
shows a rate of 30 m3/day for the Ist cycle and 18.5 m3/day for the 2nd cycle. These are 
considerably higher than the average Cold Lake operations (10.0-12.0 m3/day) reported 
elsewhere (Farouq Ali, 1994). For high early cycle oil production, the CDOR fell below 
the average method value in the 5th cycle. The effect of the oil viscosity difference is 
greatest in the first two cycles. Consequently, much less water was produced by the direct 
method during these two cycles. The WOR jumped from 2.5 in the 2nd cycle to 5.5 in the 
3rd cycle. A lower water viscosity computed by the direct method aggravated the situation 
further which is seen in Figure 6.31. 


The results and analyses presented above clearly show that the conventional direct 
method of oil/water viscosity calculation produces inaccurate and misleading results. 
Because of the extremely important role played by the oil viscosity in flow calculation, 


utmost care should be taken to estimate its value correctly. 


6.5 Conduction Between the Two Zones 


In radial analytical models, most authors (Boberg and Lantz, 1966; Jones, 1977; 
Gontijo and Aziz, 1984; Sylvester and Chen, 1988) used a single heated zone around the 
wellbore. The temperature in the heated zone at the end of the injection period was assumed 
to be the saturated steam temperature. The temperature beyond the heated zone was 
assumed to be the initial reservoir temperature. This problem was then solved in both the 
vertical and the radial directions as a hot body initially at a constant temperature cooling in 
an infinite medium through conductive heat transfer only. In warmer reservoirs, the step 
function temperature profile may not contribute any error beyond the uncertainty of the 
reservoir input data. For Canadian heavy oil reservoirs, where the initial formation 
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Figure 6.27—Effect of different methods of oil viscosity 
calculation on cumulative oil production. 
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Figure 6.28—Effect of calculating oil viscosity without areal average 


temperature on oil produced at each cycle. 
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Figure 6.29—Effect of calculating oil viscosity without 
areal average on cycle oil-steam ratios. 
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Figure 6.30—Effect of calculating oil viscosity without areal 
average on cycle calendar day oil rates. 
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Figure 6.31—Effect of calculating oil viscosity without 
areal average on cycle water-oil ratios. 
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temperature varies between 13°C and 18°C and the typical steam injection temperature is 
about 300°C, a step function temperature profile at the hot zone boundary would create 
considerable error in the heat loss calculation. One way to overcome this problem is to 
create a warm buffer zone between the hot zone and the cold outside zone. Seba and Perry 
(1969) took a similar approach for a radial model and considered average temperature 
changes in a smaller cold zone as opposed to an infinite radius. They maintained the hot 
zone temperature at the saturated steam temperature for the entire project time. Obviously, 
this would predict a highly optimistic oil production for most cases. Arthur et al. (1991) 
took the right approach in considering a hot and a warm zone in their heat calculation of an 
elliptical drainage system. The rate of heat loss from the inner hot zone during soak is very 
high. As there is no fluid flow during soak, the only mode of heat transfer is conduction. 
The contribution of conduction during the early stage of production between the two zones 
is also quite significant. Arthur et al. (1991) ignored this important aspect of the heat 
balance and considered only convective heat transfer between the hot zone and the warm 
zone. Heat loss to the over- and underburden was through conduction only. A new set of 
equations including conduction between the two zones was developed (Appendix D) in the 


current study. The importance of the conduction term is presented in this section. 


Figure 6.32 shows a comparison of the average hot zone temperature profiles of three 
cycles during production. It can be seen that the temperature at the beginning of production 
is higher for no conduction heat loss to the warm zone at each cycle. The curves are flatter 
too. On average, the difference in temperatures for the two cases through the entire 
production period is more than 20°C. Due to the higher temperatures calculated, the model 
with no heat conduction between the two zones produced 27% more oil (Fig. 6.33). Most 
of this oil was produced in the first two cycles when the effective conductive heat transfer 
between the two zones was at a maximum. Figure 6.34 shows a sharp drop in cycle oil 
production after the 2nd cycle. The cumulative water production was also 27% higher in 
the no conduction model after 673 days of production (Fig. 6.35). 


Favourable OSRs (Fig. 6.36) and CDORs (Fig. 6.37) are predicted as a consequence 
of the higher average temperature in the hot zone. The effect of excluding the conduction 
term on WORs is not as acute. The effect of high oil production in the first two cycles kept 
the water production somewhat lower and a moderate jump in WOR can be observed after 
the 2nd cycle. Figure 6.38 shows the changes in WORs for both cases. 
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Figure 6.32— Effect of conduction heat loss on average 
hot zone temperature. 
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Figure 6.33—Effect of conduction heat loss from the hot zone 
to the warm zone on cumulative oil production. 
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Figure 6.34—Effect of conduction heat loss from the hot zone 


to the warm zone on cycle oil production. 
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Figure 6.35—Effect of conduction heat loss from the hot zone to 
the warm zone on cumulative water production. 
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Figure 6.36—Effect of conduction heat loss from the hot zone 
to the warm zone on cycle oil-steam ratios. 
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Figure 6.37—Effect of conduction heat loss from the hot zone to 
the warm zone on cycle calendar day oil rates. 
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Figure 6.38—Effect of conduction heat loss from the hot zone 
to the warm zone on cycle water-oil ratios. 
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Tables 6.4 and 6.5 show the changes in the hot zone average temperature during the 
soak period with and without conduction heat loss, respectively. In Table 6.5 we can see 
that the hot zone average temperature dropped about 5 to 7°C during the 15 day soak period 
in all cycles. The heat loss was only to the over- and underburden through conduction. The 
warm zone, which is a large area, did not gain any heat from the hot zone during this 
period and the average temperature dropped from 1 to 4°C due to heat loss to the over- and 
underburden. When the conduction term is included in the heat balance, there is a 
significant temperature drop in the hot zone during soak. Table 6.4 shows a 21°C average 
temperature drop in the 1st cycle and a drop of 12°C in the 5th cycle. Although the initial 
temperature before soak is around 190°C, the heat transfer rate decreases with increasing 
cycle number due to the inclusion of the total project time in the heat loss equation. This has 
been discussed in Chapter IV. The warm zone average temperature increased by 6°C in the 
Ist cycle and 2.5°C in the Sth cycle. Considering the large volume of the warm zone, this 
indicates significant heat gain during the soak period. This temperature increase did not 
enhance the flow rate from the warm zone, because the highest temperature is only 33.8°C, 
but heat losses from the hot zone reduced the flow rate considerably from that zone. 


6.5.1 Relative Error Compared to Erroneous Viscosity Calculation 


The exclusion of conductive heat transfer between the hot and the cold zone produced 
27% more cumulative oil after 673 days of production. The effect of the “Direct” method of 
calculating average viscosity on oil production is more sever. The lower values of oil 
viscosity, calculated by this method, produced about 46% more total oil than predicted by 
the BP simulation during the same production time. This shows relatively the more 
importance of proper evaluation technique of viscosity; which is, in fact, the most 


important parameter in thermal recovery flow calculations. 


6.6 Sensitivity Study 


Although the new model successfully matched the BP numerical simulation results, a 


sensitivity analysis has been made to determine the limitations of the model. Both physical 


83 


; 1. a 
ait: nodes oranges - Art a 
on oom tae nails aaa oe 
teres dw gab Re aaa" oni a 
ssid NNER che cverers ath ieee dee sean at ‘vine = sont asd 9 
Lids ynivh anne sok a oni teelt glee dany so tab al ta 
het, wept oe ot egal detail Gt mo Yewt soaoreh eeagysrits rong 4 oes 
gue 2cled teed 94 me oye Pe & rive ¢ ioe aed ants vt 
sewn OPS 6 Genie ahuT alge Ra Re ae set 1 portly id sqm nO 
ey ft) ageok’. Se, a! Mp oy A iuineal nto get pike provi ; 
itil dw ented ake ene eit SO Beast iia 
| ( bins 2 cea ee i svie TER tect gl We povevtnt ar 

isa D6 qa Pateeige eae pee ee ore Gree Land | ees) 
crap alt le coaplel nyt areas 4 sega ut am wt ne i 


> 


inhi Silanes F porns? Sean a sik gi. hort ves" siting 
Tf ar) We es "y ey | iw) * , aS - a . Pa ce od wei at] ac Ww af orl s i 
“—s 
us ant j (Hidaka) oleh out Sees eee nish sa chron wazol 


om” enol’ » ero oF tbererimo “a aka 

¢ ek Shy wt ite ulin Gee ah eto 

.T “re = rh wy Oise wTe ate Asing so ttle BRO”: 
wih ill 26% oleae sian te aoeeee 
carte inks rege eS sions tearbeng 
, Lh it ape: er eee sa tvouahoamy 2 a oe 
by ic ao tiple gnc a aapleiod 


Table 6.4 —Average temperature change (°C) during soak with conduction. 


Cycle 1 2 3 4 5 
No. 
Hot Soak (Beg.) 198.31 191.89 196.02 190.45 188.59 
Zone 
Soak (End) 177.67 17532 180.31 LW EPs: PIG. 2 


Warm Soak (Beg.) 27.93 26.92 DM Tle) 26.50 25.91 
Zone 
Soak (End) 33.88 31.05 31.05 29.37 28.46 


Table 6.5—Average temperature change (°C) during soak without conduction. 


Cycle 1 2 6) 4 = 
No. 
Hot Soak (Beg.) 198.31 189.95 192.6 185.66 202.26 
Zone 
Soak (End) 191.98 184.1 186.85 180.25 196.45 


Warm Soak (Beg.) 27.93 26.21 26.96 2D.) ZouLs 


Zone 
Soak (End) 2/207 25.40 26.12 oan 24.71 
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and operational parameters are selected for examination. The first five cycles of the adjusted 
LGTB model are used as a base case to measure sensitivity against. 


6.6.1 Hot Zone Boundary Temperature 


The hot zone boundary temperature was selected to be 100°C. Pethrick et al. (1988) 
argued that at that temperature, heavy oil in the Marguerite Lake had a satisfactorily high 
mobility. The choice of this temperature is quite arbitrary. One has to look at the relative 
permeability, absolute permeability of the reservoir and oil viscosity-temperature 
relationship to decide upon a boundary temperature for a given reservoir. Even then one 
can argue about the validity of such a decision. 


Accordingly, the base run hot zone boundary temperature was fixed at 100°C. Four 
different temperatures — 80°C, 90°C, 110°C and 120°C — were studied to see the effect of 
the variation of temperature. Figure 6.39 depicts the cycle oil production for each case. 
Lower temperature settings (80°C, 90°C) produced much less oil in the early cycles than the 
base case, but picked up production in the later cycles. This is expected as the hot zone 
gradually gets heated with increasing cycle number. On the other hand, higher temperature 
settings (110°C, 120°C) consistently produced more oil in each cycle, with decreasing 
production at later cycles. In the case of 110°C there was an unexpected jump in production 
in the 4th cycle. This could be the result of a combination of saturation and average 
temperature calculation at the end of the 3rd cycle. The overall project performance can be 
seen from the cumulative oil production diagram (Figure 6.40). Total oil production from 
the 80°C setting was considerably lower than the base run. A setting of 90°C produced 8% 
less oil than the base case; the 110°C setting produced 9% more oil (748 m3) than the base 
run. A further increase of the hot zone boundary temperature (120°C) produced only a 
nominal increase in cumulative oil (1%) at the end of 5 cycles. From the general trends of 
the curves it can be said that any temperature setting between 90°C and 120°C of the hot 
zone boundary would give very similar ultimate recoveries. The difference in the 
performances would be in the early cycle production. A temperature of 100°C seems to be a 
balanced choice for this particular case. Figure 6.41 shows the water produced at each 
cycle. A steady increase of water production can be seen with increase of the boundary 
temperature. Incremental total production of 3,000 m3 can be seen with each increase of 
temperature up to 110°C. Additional production is about 1,800 m? for the 120°C setting. 


These results are shown in Figure 6.42. 
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Figure 6.39—Oil produced at each cycle for various hot 
zone boundary temperature values. 
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Figure 6.40—Cumulative oil production for various hot 
zone boundary values. 
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Figure 6.41—Water production at each cycle for different hot 
zone boundary values. 
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Figure 6.42—Cumulative water production for various hot 
zone boundary temperature values. 
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Without going into individual analysis of the performance parameters for each case, the 
general trends are discussed here. The lower temperature cases have better OSRs in the 
later cycles. The OSRs in the first two cycles increase steadily with increasing boundary 
temperature (Fig. 6.43). A similar trend can also be seen in the CDORs (Fig. 6.44). Water 
cuts are much higher in the first two cycles for the lower temperature settings. Figure 6.45 
shows the WORs. It can be concluded that the benefits in the early cycles for a higher 
temperature setting are traded off with lower efficiency in the later cycles. With an 8 to 9 
percent difference in total oil production at the end of 5 cycles, the sensitivity of the model 
to the setting of hot zone boundary temperature is very little, provided the temperature is 
not set lower than 100°C. 


6.6.2 Soak Time 


A soak period is unavoidable in the field because of the time needed to change the well 
from the injection phase to the production phase. Different authors have expressed 
contradictory opinions about the benefit of a soak period. Chen (1987) suggested a 
minimum possible soak period for the best oil production. Niko and Troost (1971) 
experimentally found that the oil production was insensitive to soak time. Trebolle, De Paz 
and Martinez (1993) also found that soak time had very little influence on recovery. With 
the help of a three-dimensional numerical simulation they varied the soak time from 15 to 
360 days and concluded that an inactive well did not compensate for the deferred 
production. Ferrer and Faroug Ali (1977) observed that soaking produced more oil for a 
high viscosity oil compared to the no-soak condition, whereas the situation was reversed 
for a lower viscosity oil. Resistance to the flow of oil because of water and steam saturation 
around the wellbore outweighed the benefit of viscosity reduction in the latter case. On the 
other hand, field data and numerical simulations have shown the benefit of an optimum 
soak period. Optimum soak period theory is based on the assumption that all of the steam 
does not condense right after the injection period. A few days are needed for the steam to 
condense and transfer its latent heat to the reservoir. On the other hand, this period cannot 
be too long. Otherwise too much heat would be lost to the over- and underburden. Adams 


and Khan (1969) reported an optimum period of 9 days in a California reservoir. 


The present model assumes immediate condensation of the steam after the injection 


period and does not account for any steam distillation. In a cold heavy oil reservoir, this 
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Figure 6.43—Oil-steam ratio at each cycle for various 
hot zone boundary values. 
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Figure 6.44—Calendar day oil production rate at each cycle 
for various hot zone boundary values. 
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Figure 6.45—Water-oil ratio at each cycle for various hot 
zone boundary values. 
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would not be too unrealistic. Different soak periods between 0 to 18 days were considered. 
As shown in Figure 6.46, cycle oil production steadily decreased with an increase of soak 
time up to the 2nd cycle. After the 3rd cycle, higher soak periods (15, 18 days) produced a 
little more oil than the lower soak times (5, 9 days). This is due to more available oil and 
late heating of the formation. A better picture of the performances can be seen in Figure 
6.47 where cumulative oil productions are displayed. A no-soak time produced the highest 
amount of oil at the end of 5 cycles. This is expected with the assumptions made for the 
model. Total oil production always declined with the increase of soak time but for a wide 
variation of soak period between 5 to 18 days the decrease was only about 5%. Except for 
the no-soak line, all the other curves in both cycle (Fig. 6.48) and cumulative (Fig. 6.49) 
water production diagrams clustered together to see any significant difference between 
them. Although these results suggest better performance with a shorter soak time, loss of 
recovery does not increase significantly with a longer soak period beyond 5 days. To 
prevent loss of production due to shut-in, a soak period between 0 to 5 days seems to be 
the best strategy. Remembering the limitation of this model mentioned earlier, caution 
should be taken in using this model to establish soak periods in any CSS operation. 


6.6.3 Increasing Steam Volume 


Several experiments were done with the steam injection strategy. The first of these tests 
was increasing the steam injection volume. Injection periods for each cycle were kept 
constant for all the tests. The injection rate was varied from 150 m3/day (cold water 
equivalent, CWE) of steam to 400 m3/day. Figure 6.50 shows a steady and uniform 
increase in oil production at each cycle. This was due to increased heat input into the 
reservoir. For each m3/day (CWE) increase in steam injection, about 5 m3 of incremental 
oil was produced in the Ist cycle. Incremental oil production per cubic meter of CWE 
decreased with an increasing amount of steam injection. An increase of injection rate from 
150 m3/day to 200 m3/day for a total of 148 days of injection time over the 5 cycle period, 
produced 32.3 m3 of incremental oil per m3/day rate increase. For the same operation 
period, the incremental oil was 27.5 m3 per m3/day for a rate increase from 300 m3/day to 
400 m3/day. Cumulative oil production for each case is shown in Fig. 6.51. The 
diminishing return with increasing steam injection rate was also observed by Martin (1967) 
and Davidson et al. (1967). In a steamflood study, Hong (1994) observed a similar trend. 
Boone et al. (1993) suggested a strategy of massive steam injection (over 2,000 m3) for 
short periods (hours) as opposed to a constant lower injection rate for longer period to 
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Figure 6.46—Oil production at each cycle for different soak periods. 
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Figure 6.47—Cumulative oil production for different soak periods. 


92 


. ; = > i 
WwW 7 ieay 
A o ss 


shu 
ion SEs Tog orth’ 1CF GinED (see er Stteageeany DSB 
; ae | 70 


1.240% 


| —o— 0 days 
} — -5days 
} — @ -9days 
1 104 - -* -- 15 days 
} - -+- - 18 days 


8000 


6000 


4000 


Cycle Water Produced, m® 


2000 


Cycle No. 


Figure 6.48—Water production at each cycle for different soak periods. 
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Figure 6.49—Cumulative water production for different soak periods. 
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Figure 6.50— Oil production at each cycle for various steam 
injection rates. 
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Figure 6.51—Cumulative oil production for different steam 


injection rates. 
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increase the areal thermal conformance. This strategy was based on changing fracture 
orientation and geometry. Increased steam production involves higher boiler capacity and 
higher operational costs. A detailed calculation of heat loss through the tubing has to be 
included in the analysis. Economic factors play an important role in steam injection 
strategy. Decisions regarding injection rates have to be based on cost per cubic meter of oil 
produced. Cycle water (Fig. 6.52) and cumulative water production (Fig. 6.53) increased 
proportionately in a similar manner. 


The production performance parameters give an excellent perspective on the steam 
injection strategy. Whereas higher injection rates produced lower OSRs (Fig. 6.54), 
CDORs (Fig. 6.55) varied inversely with injection rate. Better OSRs were achieved 
through lower injection rates. Both 150 m3/day and 200 m3/day attained identical values, 
whereas 400 m3/day was the worst performer. Obviously any increase of injection rate 
would produce even lower OSRs. Depending on the operational break-even point, field 
Capacity and initial investment, each reservoir would have a different optimum strategy. 
Shepherd (1979) also found that increasing steam volume increased the CDOR, although 
recovery from a well was not significantly affected. Figure 6.56 shows lower WORs for 
higher injection rates. This could be explained as a result of larger contact area for 


increasing injection rates. Steam travels farther from the wellbore allowing a higher oil cut. 


6.6.4 Slug Size 


About fifty percent of the cost in thermal projects involve generating steam. It is very 
important to explore the injection strategy with a given steam generating capacity. In the BP 
simulation a total of 37, 000 m3 CWE of steam was injected with increasing slug sizes. 
Three more simulation runs were conducted with different slug sizes. Large-small-large 
(LSL), constant and decreasing slug sizes were used, keeping the total steam volume same. 


Table 6.6 shows the four steaming schemes. 


Figure 6.57 shows the cycle oil production. Both the LSL and the decreasing slug 
strategy produced the same maximum oil in the 1st cycle with slug sizes of 8,250 m3 and 
10,000 m3 respectively. With the lowest slug sizes, the decreasing slug scheme produced 
less than 300 m3 of oil for each of the last two cycles. With the decrease of slug size, the 
LSL oil production dropped abruptly in the 2nd cycle and then gradually picked up after 
that. A constant slug size produced little more oil than the increasing slug size for the first 
two cycles and then the production decreased steadily for the next three cycles. Figure 6.58 
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Table 6.6—Slug sizes (m3) for steaming schemes. 
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Figure 6.52—Water production at each cycle for different steam 
injection rates. 
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Figure 6.53—Cumulative water production for various steam 
injection rates. 
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Figure 6.54—Oil-steam ratio at each cycle for various steam 
injection rates. 
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Figure 6.55—Calendar day oil production rates for different 
steam injection rates. 
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Figure 6.56—Water-oil ratio at each cycle for different steam 
injection rates. 
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Figure 6.57—Oil production at each cycle for different slug 
sizes with constant total steam. 
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Figure 6.58—Cumulative oil produced for different slug sizes 
with same total injected steam. 
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gives a better idea of the overall performance in terms of the total oil recovery. The 
increasing slug size (BP simulation) produced the maximum amount of oil after 5 cycles, 
doing marginally better than the LSL scheme. The constant slug strategy produced 5% less 
oil than the increasing slug strategy. The decreasing slug size strategy was the worst 
performer of all the schemes with 21% less oil production. From these results, it is seen 
that with a large slug size much oil can be produced in the 1st cycle but, reduction of slug 
size reduces the subsequent cycle production drastically. It is very important to fill up the 
voidage in the reservoir with increasing slug sizes. Decreasing the slug size after an initial 
large steam volume does not expand the hot zone any more. It would fill up the already 
heated zone and just bring up the temperature without reaching new oil. Cycle water 
production (Fig. 6.59) for each of the schemes was much higher than for the increasing 
slug size scheme. Consequently all the three cases produced more water after 5 cycles (Fig. 


6.60), with the decreasing slug size strategy producing the maximum water. 


Figure 6.61 shows the OSRs for each cycle. Again, increasing the slug size shows the 
best overall performance. The same conclusion can be made from the CDORs (Fig. 6.62). 
The WORs for all four cases were fairly close except for a dismal performance again by the 
decreasing rate scheme. The WOR value rose over 30 for the last two cycles. Figure 6.63 
shows the results. With much less heat supplied, the reservoir average temperature would 
drop favouring large water production from the same heated zone created in the Ist cycle. 
Unless one is aiming for high oil production in the first couple of cycles, increasing the 
slug size is the best injection strategy in CSS operation. A close second would be the LSL 


injection with a quick pay-off in the 1st cycle without loosing much overall efficiency. 


6.6.5 Time Step 


Analytical models, from the nature of their solution procedure, should not be sensitive 
to time step sizes. In this particular model, the most important factor for flow calculation is 
oil viscosity. Viscosity for the next time step is calculated from the current average 
temperature. As the average temperature is a function of time, so is the viscosity. The 
highly non-linear relationship between temperature and viscosity may make this model 
sensitive to time step size. Gontijo (1983) reported a considerable production drop with an 
increase of the time step size using a similar viscosity calculation method. Figure 6.64 
shows the effect of different time step sizes on cycle oil production: 0.5, 1.0 and 2.0 days 


show virtually no difference. The first significant change in production pattern is seen 


101 


10! 


ts pa lychee 

gout PO Feb cholo la euinas at 3 | 

mae Sch ghar ge tiie ee ee goku att 
+1 pad. Gib anagem seat oot sgbouleng Vic tent FT i 
1g to oalvouben it sive act oie melee pes: 
1 ge DT od idiot aoe ac .Olenth eahonbeng & 
twisind ma catia seb gotta qiizeseeed arhe guile rriepsaais Hite 9 
eben aattteape, (is Hilger 0 Locate cig rine tod =) bongEe 1 
catte wlag flo, weal ge rage ‘potitw sip Ml qe ‘lek 
Leston << Wt watt outer thorn sane area, ot hg tae TS 
aslives 2 satin wiles ery deve ong tees nae wd Lip cio 
Hew Ge VME Tet ANG CLOT Bae Bae fi 


7 


¥ 
deel 
5 aD 
7 


mn 
- 


rl Ay eee fi let aes wt 2AeKD wit 249, 
1) ert ab: 1 A000 dient cat 004 eee aie cee ObT on 
ope careeemctieey Cramer fe) goa Gi yale) dale ene ae 5 
(PP dotret wad by ah eR ni het aoe OS a eee 
' + ot) Lien ett a) Readers Geel Advit av 

t dp 1: oat ¢ Geet aay) le ee) Caley ee 

re we orld Gar a) otra My git at 

Lino Levte eeapy h Alay, OS eee irae 


ee ee 2 ten 8 sate of Phe yan aad 
ry 


won Oke: Swlkedanainlte dale saan 
zi err a ee chee ea piabia aged (50h Oli iguich? 
rT Mirren ate iiinlan 0p mae) tage 
ot ort af) 4) Ge onlr eter a 
inhi ee «ely Uplaizey aia een 
on {uve See sth alee “ste ae eam MiP 
in eats daar sapee:?- REY ps 4 ese 
anon G4 bene 04 2 ie . we. ee 
seer 03 ola ihal, “a 0%) 


Cycle Water produced, m?® 


Cumulative Water Produced, m° 


3.5 104 


3 104 


2.5 104 


2 104 


1.5 10* 


1 104 


5000 


—0e— BP Simulation 

—4 _ - Large-small-large Slug Size 
— © - Constant Slug Size 

- - - - Decreasing Slug Size 


Cycle No. 


Figure 6.59—Water production at each cycle for various slug 
sizes with constant total injected steam. 
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Figure 6.60—Cumulative water production for various slug 
sizes with constant total injected steam. 
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Figure 6.61—Oil-steam ratio at each cycle for various slug 
sizes with constant total injected steam. 
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Figure 6.62—Calendar day oil production rate at each cycle for 
various slug sizes with constant total injected steam. 
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Figure 6.63—Water-oil ratio at each cycle for various slug 
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Figure 6.64— Oil produced at each cycle for various production 
time-steps. 
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when the step size is increased to 5 days. The lst cycle production increases and starts 
dropping from the 3rd cycle. An increment to 10 days still follows the 0.5 day trend 
reasonably closely. Doubling the time to 20 days changes the production pattern 
considerably. In terms of cumulative oil production, the differences for all cases are very 
small. A 20 day step produced 283 m3 less oil than the 0.5 day case after 600 days of 
production, which is less than 3.5% production loss. Figure 6.65 shows the results. Cycle 
water (Fig. 6.66) and cumulative water (Fig. 6.67) curves for all the cases are so close that 
it is difficult to distinguish one from the other. The improved response toward time step 
size is due to the adaptation of the “Average” method of oil viscosity calculation. In this 
method, the oil viscosity is not as sensitive to temperature changes as the “Direct” method. 
All the runs in the present study used a 2 day time step size. This enabled faster 
computation without losing accuracy. 


6.6.6 Steam Temperature 


Another feasible option to increase the heat input into the reservoir is to increase the 
steam temperature. The cost of producing steam at higher pressure or to superheat the 
injected steam is high. Two runs were conducted to see the benefit of high temperatures. 
Steam saturation temperature (pressure) was increased by 50°C and 70°C above the base 
run. Steam quality for each case was kept at 0.75 to see the effect of temperature. This 
would actually be disadvantageous for higher temperature, as the latent heat content 
decreases with increasing pressure. Steam enthalpy remains almost constant with less latent 
heat at a higher temperature. Figures 6.68 and 6.69 show the effect of higher temperature 
on oil production rates and average hot zone temperatures respectively for the first cycle. 
Higher temperature enabled higher rates for the high temperature cases. Figure 6.70 shows 
the cycle oil production. For the higher temperatures, the 1st cycle production was higher 
and after that the responses were somewhat erratic. Cumulative oil productions are shown 
in Figure 6.71. It can be seen that higher temperatures produced a larger amount of oil. 
Interestingly, increasing the temperature from 50°C to 70°C produced marginal incremental 
oil after 5 cycles. Even the 50°C case produced only 7% incremental oil. Although 
increasing temperature makes the oil more mobile, more heat is also lost in these cases. 
Moreover, the total heat input remains almost constant. In Figure 6.72, it can be observed 
that at a higher temperature the fluid heat loss is higher. Both cycle water (Fig. 6.73) and 
cumulative water (Fig. 6.74) follow similar trends as oil production. From these results, it 


can be concluded that increasing the steam temperature follows the law of diminishing 
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Figure 6.65—Cumulative oil production for various production 
time-step sizes. 
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Figure 6.66—Water production at each cycle for various 
production time-step sizes. 
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Figure 6.67—Cumulative water production for various 
production time-step sizes. 
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Figure 6.68—Effect of steam temperature on 1st cycle oil 
production rate. 
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Figure 6.69—Effect of steam temperature on 1st cycle 
average hot zone temperature. 
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Figure 6.70—Oil production at each cycle for different 
steam temperatures. 
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Figure 6.71—Cumulative oil produced for different steam 
injection temperatures. 
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Figure 6.72—Heat lost to fluid production at each cycle 
for varying steam injection temperatures. 
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Figure 6.73—Water production at each cycle for various 
steam injection temperatures. 
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Figure 6.74—Cumulative water produced for different steam 
injection temperatures. 
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return. Similar observations were made by other authors (Davidson et al., 1967: Gontijo 
and Aziz, 1984). Martin (1967) actually got better performance from low pressure steam 
injection. Moreover, a substantial jump in cost does not justify increasing the steam 
temperature for such nominal incremental oil. It should be noted that the aforesaid 


investigators used very different mechanistic models from the one used for the present 
study. 


6.6.7 Steam Quality 


The last controlling factors for the heat content of the injecting steam is quality. For an 
effective steam operation, sufficient quality steam at an optimum rate has to be injected. 
Better quality increases the heat content and more importantly, it increases the latent heat 
part of the steam enthalpy. Latent heat plays a significant role in steam zone expansion 
(Marx and Langenheim, 1959; Wheeler, 1969; Mandl and Volek, 1969). Steam qualities of 
0.6, 0.75 and 0.96 were used to see the effects on oil and water production. 


Figure 6.75 shows a monotonic increase in cycle oil production with increasing steam 
quality. The impact of increased heat input is greater with increased slug size. In the first 
cycle with the smallest slug size, the difference in heat input was not sufficient to change 
the average hot zone temperature. Cumulative oil production also increased monotonically 
with increasing steam quality. Figure 6.76 shows a 10% increase in total oil production for 
an increase in quality from 0.6 to 0.96. This improvement is not as much as in California 
(Gontijo and Aziz, 1984) or in Venezuela (Trebolle et al., 1993). This is due to the high 
injection pressure in the Cold Lake area (over 10 MPa). As the pressure increases, the 
contribution from latent heat decreases and the effect of improved quality diminishes. Cycle 
water production (Fig 6.77) slightly increased after the 2nd cycle for the increasing steam 
quality cases. There was no significant difference in the total water production (Fig. 6.78). 


All three production performance parameters — OSR (Fig. 6.79), CDOR (Fig. 6.80) 
and WOR (Fig. 6.81) show favourable results for higher steam quality. This could be 
somewhat misleading, as none of these parameters considers the cost of improving the 
steam quality. Performance has to be assessed on the basis of economic factors involved in 
steam production. Wellhead steam qualities of 0.75 to 0.8 are most widely produced in the 
field. 
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Figure 6.75—Oil production at each cycle for various 
steam qualities. 
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Figure 6.76—Cumulative oil production for various steam 
qualities. 
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Figure 6.77—Water production at each cycle for various 
steam qualities. 
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Figure 6.78—Cumulative water production for various steam 
qualities. 
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Figure 6.79—Oil-steam ratios at each cycle for various steam 


qualities. 
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Figure 6.80—Calendar day oil rate at each cycle for various 


steam qualities. 
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Figure 6.81—Water-oil ratios at each cycle for various steam 
qualities. 
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6.6.8 Initial Reservoir Temperature 


To a considerable extent, the effectiveness of steam processes depends on the initial 
reservoir temperature. While the average reservoir temperature in California varies between 
30 to 40°C, Alberta heavy oil fields have formation temperatures between 10 to 20°C. For 
many Venezuelan reservoirs this value exceeds 50°C. In the present study, the same 
viscosity-temperature relationship was used for all the runs. This relationship actually 


varies from reservoir to reservoir and should be considered for a true comparison between 
reservoirs. 


Initial reservoir temperature determines the amount of net heat input into the reservoir. 
With a high initial temperature, the available enthalpy would be less as the sensible heat 
transfer from the steam decreases. The effectiveness of steam injection varies directly with 
the available enthalpy to the reservoir. Figure 6.82 shows the net heat injected into the 
reservoir for different initial reservoir temperatures. As a consequence of greater heat input, 
lower initial temperatures would have a higher hot zone average temperature. This would 
allow more production. Figure 6.83 shows the effect of initial reservoir temperature on 
cycle oil production. Up to the 20°C case the cycle oil production response shows little 
change. Faster growth of the hot zone in the early cycles for higher initial reservoir 
temperatures (40°C and 60°C) produced saturation values quite different than the lower 
temperature cases (10°, 16° and 20°C). Because of these saturation values and lower 
production rates, the cycle oil production curves for 40°C and 60°C follow different trends 
than the lower temperature cases. Figure 6.84 shows a reduction in cumulative oil 
production when the initial reservoir temperature is increased beyond 40°C. For the 10°, 16° 
and 20°C cases the values are almost identical. One advantage of a higher initial temperature 
is quicker early growth of the hot zone. Figure 6.85 shows the progression of the 100°C 
isotherm. At Tp=60°C, the 100°C isotherm moves 0.6 m more beyond the fracture face 
after 5 cycles. This gain was not sufficient to offset the effect of a lower average zone 
temperature. Water production for the higher temperature cases (40°C and 60°C) was also 
less due to the lower average temperature. Figure 6.86 displays the cumulative water 
production up to the 5th cycle. It is clear from these results that due to more available 


enthalpy, steam injection is more effective for colder reservoirs. 
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Effect of initial reservoir temperature on net 


heat injection. 
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Figure 6.83—Oil production at each cycle for different initial 
reservoir temperatures. 
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Figure 6.84—Effect of initial reservoir temperature on total 
oil production. 
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Figure 6.85—Effect of initial reservoir temperature on hot 
zone propagation. 
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Figure 6.86—Effect of initial reservoir temperature on total 
water production. 
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6.6.9 Formation Thickness 


Oil recovery from steam stimulation is affected by the pay thickness. The effect would 
be more prominent due to gravity drainage. The capability of the model to predict this effect 
was investigated. Reservoir height for the base case was 23.8 m. Two more runs with pay 
thicknesses of 30 m and 45 m were simulated. 


Cycle oil productions are presented in Figure 6.87. As expected, a greater pay 
thickness produced a larger amount of oil. Interestingly, increasing the thickness from 30 
to 45 m did not produce any significant amount of incremental oil over the 5 cycle period. 
Results are shown in Figure 6.88. This can be explained due to the fact that the injected 
steam volume was kept constant for all three cases. For the initial increase in thickness a 
reduction in injection rate per meter of pay was outweighed by the lower heat loss. 
Although heat loss to the over- and underburden reduced further with the increase of 
thickness, the same amount of heat was not sufficient to heat up the entire formation to a 
higher temperature level. As a result the oil production did not improve. Heat losses to the 
adjacent strata for the 45 m thickness was 45% to 50% less than for 23.8 m. Heat losses 
for all three cases are shown in Figure 6.89. Apart from a lower heat losses in the thicker 
formations, gravity flow was also more effective. Water production for thicker formations 
was much higher too. Results are shown in Figures 6.90 and 6.91. 
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Figure 6.87—Oil produced at each cycle for different pay 
thicknesses. 
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Figure 6.88—Cumulative oil production for different pay 
thicknesses. 
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Figure 6.89—Effect of formation thickness on heat loss to 
over- and underburden. 
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Figure 6.90—Water produced at each cycle for different pay 
thicknesses. 
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Figure 6.91—Cumulative water production for different pay 
thicknesses. 
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Chapter VII 


CONCLUSIONS 


An analytical predictive model was developed with the objective of defining heat and 
fluid flow in heavy oil CSS operations more accurately. Contributions made by this study 
and the conclusions reached are as follows: 


1. For the first time a comprehensive analytical study of the elliptical shaped steam 
zone with a wellbore at the centre has been made. Flow into a circular 
wellbore in the centre of an elliptical flow geometry does not represent the flow 
pattern in the field (heavy oil). It underestimates production. 


2. This is the first analytical model which considers gravity flow into a partially open 
fracture. It has been compared with a simulation study and yields satisfactory 


results. 


3. An area weighted true average oil and water viscosity evaluation technique has been 
developed for the first time. It has been shown that the conventional method of 


calculating viscosity may produce rates as five times too high. 


4. The importance of horizontal heat conduction in a two-zone heat balance has been 
established. Inclusion of the conduction term provides improved average hot and 


warm zone temperatures. 


5. The model shows plausible response to variations of a wide variety of physical 
and operational parameters: steam volume, slug size, steam temperature, steam 
quality, formation thickness, hot zone boundary temperature, initial reservoir 


temperature, soak time and production time step. 


6. This study provides a quick means of assessing CSS operations and optimization 


and providing a sensitivity study. 
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7.1 Suggestions for Future Studies 


Several ideas were explored for further studies in the process of developing this model. 
The more important ones are: 


1. To investigate the effect of a horizontal fracture. 


2. To accommodate the concept of dilation and compaction, especially in the closed 
fracture case. 


3. To implement the average temperature calculation method in a horizontal well steam 
injection process. 


4. To apply the model to more field data. 
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Appendix A 


CORRECTION OF SYLVESTER AND CHEN (1988) EQUATION 


Sylvester and Chen (1988) developed a steam soak model for a homogeneous gravity 
drained pressure-depleted reservoir. They argued that the proposed gravity drainage model 
by Towson and Boberg (1967) had several uncertain terms. One has to calculate fluid 
heights at the wellbore, heated radius and the average fluid height in the heated zone. Chen 
(1987) in his Ph.D. dissertation derived an equation which eliminated all of these 
parameters and was replaced with one parameter. This flow equation was derived based on 
Equation A7 of Towson and Boberg (1967) which is Equation III-19 in Chen's thesis. 
There is a serious error in the derivation of this equation. Chen's Equation III-19 is 
Equation A.1 below. Starting with this equation, Chen's error will be shown: 


Oe ae eae 
eesti hp —h? (Aan) non 
(R=? yn S a5) 


w 


If B= fin _ a5 then Equation A.1 can be written as 
. 


w 


a 42 
hah. aS (A.2) 
After rearranging, 
h? —-8hBh, + 8h B-h? =0. (A.3) 


This is a quadratic equation in hp, so the solution is 


Fe T, p2 7,2 2 
a ShB + 64h B? -32h B+4nh? Ces 


h 2 Y 


h, = 4hB+\16h B? + h2 -8h' B. (A.5) 


138 


Chridy AIRSET 4 na ecw aap as 
Pe . rt ee ‘iy, LRP 5 =v bial. is pha zy 
4 egriee. oe acinomae pee 
‘“” { ‘ | pi 7 oe M wt aqeetiaads tena adhe : Pra - 
we ek 
} pref ayia’ wo Mate ah + - ' 
RD Hy Sp eagy Hrmmaec! 0% afalel 

ee. ots _ te wv Een) eats 
7 abt, 4 Le tf. a (ere bi aedetl ba o), 2 ' 

in a tial rt naiimeriab oipnt te ante 
rete Hie ond Got? 2288 rors a 


jt) ee 


| 
ohne 


A? 


BAY 


qh) 


Instead of the (+) sign in front of hy2, Chen proceeded with a (-) sign and derived a 


relationship between hy and hy in terms of B which is (using his relationship 
h=(h, +h,)/2): 


2B+1 
h = ie 


If one proceeds with the right sign, the above relationship cannot be reached. Taking 
squares of both sides of Equation A.5 and assuming that h =(h, +h,)/2, as Chen did, it 


can be written 


2 2 
W2 — he = a, (AS) ¥ (Hethek | (A.7) 


hy, — hy, = 2B(hy, - hi). (A.8) 


This gives us B=0.5 or with rearrangement hw = hh, neither of which gives us any 
meaningful correlation. Chen (1993) was apprised of this error and recognized the mistake. 
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Appendix B 


FLOW INTO A WELL LOCATED AT THE CENTRE OF AN 
ELLIPTICAL DRAINAGE BOUNDARY 


The standard procedure to solve this problem is to employ elliptical coordinates. Van 
der Ploeg, Kirkham and Boast (1971) presented a closed-form solution for steady-state 
saturated flow into a finite circular wellbore in an elliptical confined aquifer. They solved 


the flow problem for a number of well locations in an isotropic and homogeneous aquifer. 


B.1 Well at the Centre of an Ellipse 


The equation of ellipse in a rectangular coordinate is 
Aes). 
Si (B.1) 
a b? 


Using polar coordinates (R,®), x and y may be written as 
x =Rcos0, (B.2) 
y =R sinO, (B.3) 


and the equation of ellipse for a point P(R,®) at the boundary is 


2 «2 
re 2 Oat. (B.4) 


B.1.1 Boundary Conditions 


a. D=0 forr=Iw 0<O0<7n/2 (B.5) 
b, D= 1 Tout = hs 0<O0<7n/2 (B.6) 
c. 0/00 =0 for 9 =0 Iw <r<a (B.7) 


d.d@/00=0 forO0=n/2 tw<r<b (B.8) 
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Laplace's equation in polar coordinates is 


CHU Ae B 
pcg | pas lal el |; (B.9) 


The solution of this problem should give an expression of @ which would satisfy the 
boundary conditions and Laplace's equation. 


B.1.2 Solution 


The authors used the Gram-Schimdt method as modified by Powers, Kirkham and 
Snowden (1967) to determine the solution as 


N 
D= PV Ay lin(19), (B.10) 
m=0 


where 
= Oe 2. oN: 


N= 0; 102; .. 00: 
ate 2m 
a (ar 
“a im 2m 
a 


Replacing r with R at the boundary and using Equation B.4 to express R in term of 0, 


be (r,.0) = cos2mé. (B.11) 


u(r,6) may be written as 


(F) [ee] 
=a |) ead pcarsey 
We Oye eee Le TS cosa), (B.12) 


where 


2 2 0 = 
pe -( ore se (B.13) 
a b 


TY (¢) rn ae Gy 
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The hydraulic head may be written as 


N 
De te 2, Arn'tn(8). 0<O<n/2 (B.14) 


Applying boundary condition (b) the hydraulic head may be written as 


N 
Dae=ii-a) Aue < 1056 < 1/2 (B.15) 


m=0 


Powers et al. (1967) derived a table of orthonormal functions to solve potential flow 
problems like seepage of steady rain through soil bedding. In accordance with Powers et 
al., the two constants to determine An» in this problem are 


Wo 


(1)u,,(0)d0, m=0,1,2,...,N (B.16) 


on la 


3 
Uny = | Up (Ou, (O)d0. m=0,1,2,...,N ns<m (B.17) 
0 


The parameters u,,(@) and u,(@) may be determined from Equation B.12. When N 


—>ee, boundary condition (b) is satisfied exactly. It is difficult to calculate the above 
integrations analytically. A more accurate Gaussian quadrature numerical integration 
method was used to evaluate the functions Wm and Umn, instead of the Simpson's rule 
suggested by the authors. It can be seen that terms with zero subscripts produce 
indeterminants. Using l'H6pital's rule ug would be 


in 


Uy(8)=— (B.18) 
In— 
x 


w 


Therefore, the hydraulic head may be written as 
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After the values of Wm and umn are determined, all An values can be calculated by 
using Table 2 of Powers et al. (1967). 


cos2mé@. (B.19) 


B.2 Sequential Formulas Presented by Powers et al. 


‘ao Oe ee Seve (B.20) 
jane 

Pea GD iO 1s 2eN (B.21) 

ip = Sees m= 0, 1,27.,N (B.22) 
n=0 

lig = 16 Sau, Teale ot ss (B.23) 

oct 
ie Oem Sous i= 203, 480 (B.24) 


oe UR AG eet Bg 


G. ae m = 1, 2, 3,4,... (B.25) 
T= Oo 12 3 pakey Oi t 
m-1 
Deeetine-3)C,, Dy m= 0, 1, 2, 3, 4, ... (B.26) 
n=0 


= Onl e254 
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B.3 Well Discharge 


Using Darcy’s equation for potential flow for unit thickness, the flow rate q can be 
written as 


Ww 


@ 
a=-Kao (2 ) Zieh oe, Sti), (B.27) 


r=K, 


where K is the hydraulic conductivity. Multiplying the above equation with the thickness h, 


the total flow is 


=-Knae|(2) rd@ reso) (B.28) 


r=n, 


The flow rate for a well in the centre of an elliptical drainage system is 


2TKhAy AP (B.29) 


a 
In— 
: 
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Appendix C 


HEAT LOSS CALCULATION TO THE CAP ROCK 


Vinsome and Westerveld (1980) proposed a simple solution to the thermal diffusivity 
equation in the cap and base rock. A fitting function 


Tice (7, + pz+ Tile (Coly 
was used to approximate the temperature profile in the cap rock, where To is the interface 
temperature and d is the diffusivity length. The parameters p and q are to be determined. 
This is a sensible function to pick to describe the temperature profile. The function 
automatically satisfies the two boundary conditions (T=0 at z=0 and T=0 at z=-). 


The equation of heat flow into the cap rock is 


or OT 
eam oo ee 


The above equation can be written satisfying the interface condition as 


2 
z=0 


Inserting Eqn. C.1 into Eqn. C.3, together with a finite-difference discretization they 
came up with a solution for the heat loss rate and the energy stored in the cap rock. These 
equations were presented in the main text (Eqns. 4.16, 4.17) and are repeated here for 
convenience. 


ii 
d= (2-9) 


E. = “d(T, + pd +2qd’) 


h 
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Arthur et al. (1991) defined the parameters in the following manner, 


Ty At ity ae a\(T, = T") 
d a, At 


Ne 3d° + a,At 


d*(T, -T’) 


2pd—-T, + 
p : a, At 


and 
I" =T,'d" + p"(d")’ +29"(a"y’ 
The heat loss term was expressed as 
Gy = Ty ~ kOe, 
where a, = i( = a, and 


a, At d° 
LEIS) Se oa ey a ree n=l 
d(3d? +a,At) @,At(3d? + a,Ar) 
(2 ees ee ae 
" @,At(3d? +@,At) 3d? +0,At" 


(C.4) 


(C.5) 


(Ee) 


(Cx7) 


(C.8) 


(Cy) 


(C.10) 
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Appendix D 


AVERAGE ZONE TEMPERATURES DURING SOAK AND 
PRODUCTION 


A heat balance is done on two concentric elliptical cylinders (Fig. D.1) to obtain the average 
temperatures. The parameters A; and A? are the heat conducting areas of the hot and warm 
zones adjacent to the cap and base rocks. They are assumed to be same at the top and at the 
bottom. The parameter Ap is the area between the two zones. 


Figure D.1—Physical diagram for the heat balance. 


Heat balance in the hot region (zone 1) 
7) 
~ 5, (AAMT) 7 GAC A(T; x T,) i IPrC p(T, > T;) +0,(1)+O,(). (D.1) 
Heat balance in the warm region (zone 2) 
7) 
~ 5, (AshMT,) = DePaC pati -T,)+Q,(t)-Q,(¢). (D.2) 


Assuming the flow rates and the fluid and rock properties to be constant during a time 


step, the above equations can be written in expanded form as 
fy C.(T, -T,)+(2A,+4 Tak i3 
~A,hM, = GPC ll —Tr)— 92P2 p2(Ty — ,) + (2A, + » (OT; — KO%,,), (D.3) 


dT. 
—A,hM, A = GP ,Cyo(T, — Ty) + 2A, (OT, — Klay) -— Ay (OT, - kO,;). (D.4) 
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The hot zone equation can further be expanded into 


alee DPC pr (Gone IP2rC a(2A, + A,) T + k,y,(2A, + A,) 


T. - xP ’ 
dt AhM, AhM, UD, AhM, A hM, 


Ml __5 (1, -1,)+ By(t-T,) B+ Be 

A =(B, ~ B - Bs) + BT_— Bil, + B, 

qd, _ pr 

ee RE BT, + BT + Bs: 

1 = BI, Bis + By. eed 


Using the operator D for d/dt, the equation may be written as 


(D — B;)T, a sate =p, =U he 
Similarly for zone 2, 

aT, 

mae ee — B,)T, —(B, =fiig)be + B, — By 

a = Biot, ai B37, ce ve a 


Using the operator D 


BysT, + (D - By )T, — By, = 9- (D.8) 


The constants are defined as 


B= DPC (D.9) 
' AAM, ” 

ee MPrCp2 é (D.10) 
> A|AM, 

ae a.(2A, + A,) (D.11 
3 = ’ 


AM, 


: ea 7 7 
i ae i. s f 
aan 7 


. 
4 att Baten deta OF 24 Sahl iat 


wall beens 9 sa i= 


° mm! 
a. 


i: oy ae 
(“whe 
' 


: . * oe 
io sd Gen oC apr oN oe See 


2 k,0t(2A, + A,) 


Ps ALM 

B; = B, - B, - B;, 

Bs = BT, + B,, 

5 es 
A,hM, 
2. 

B, ~ iM,’ 
ZK. 

eevee 

Bi = B, - B, 
A,Q. 

Bi ae 

fi, = an, 
A,hM, 

Bs = B, - Bu, 

Bis = By — Bir- 


The following differential equations have to be solved simultaneously. 
(D- B;)T, + BT, = By 
BysT, +(D- Byo)T = Bs 
The determinant of the matrix 


(D-B;) B, 
B,3 (D-B,,) 


is D2-DA-B, where, A=B5+f10 and B=B13B2—-B1085. Now it can be written 


(D.12) 


(D.13) 


(D.14) 


(D.15) 


(D.16) 


(D.17) 


(D.18) 


(D.19) 


(D.20) 


(D.21) 


(D.22) 


(D.23) 


(D.24) 


(D.25) 
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B, B, 
AT, = = DB, -P, 
Bi, (D- By) P (D.26) 
and 
(D— B,) B, 
AT, a =) 14 al as : 
Be [bon P D 2, 


where, P=B6B10+B2B14 and Q=B5814+B6813. The characteristic equations for Ty and T2 
are found to be, 


if = Ce" Se C,e""" +, (D.28) 
= At At OQ 
Ty =Goet + Ce ey (D.29) 
where, 
A+7VA’+4B 
A ri (D.30) 
A-VA?+4B 
Ly eee (D.31) 
2 
A=; + Bio; (D.32) 
B= B,,B, — BsB,o. (D.33) 
P= BsBio z BBs (D.34) 
and 
Q= BS B,4 a Be Bis: (D.35) 


After applying the initial condition and simplifying, the two average temperatures can 


be expressed as 


T, = Ce" + Ce" + 7 fobs) 
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Co DKS) 
B, B, awe ee 
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0 E 
B,T, =k (B, —1,)- B, OAs 

is , (D.38) 

A, ion A, 
CC. eek (D.39) 


Appendix E 


SATURATIONS 
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The mobile oil/water saturations (Arthur et al/., 1991) at.any time n+1 in the hot and 


warm zones are given by 


cer = ge 


omh.w 


1 


ie — cn 


wmh,w 


1 
where for hot zone, 


oi aa? Qi i Qo2 1 Qw2 


xX, = - 
ae 
— ~Fo1Po1 F Fo2P 02 

ee ee 
PoY hot 

xX,= —QyiPwi a3 GwoPw2 

3 ae . 

PasVivs 


For the warm zone, these variables are 


be x Qo1 ay = Qo2 = Qw2 é 


me ; 


warm 


Saturations are computed from the following equations 


Nae) = Se ihe os Does 
Ne ee — 1 ie inet oe 


omh,w exp ae sa ae “(1 = exp | ee Si ) ), 


& n+l _4n Bys wes ntl _ yn 
amhwe oe " : 1+ 3(1—exp a | Ih 


(E.1) 


(Ee) 


(E.3) 


(E.4) 


(E.5) 


(E.6) 


(E27) 


(E.8) 


(E.9) 


(E.10) 


mt 


Jo md Aad On awerte OO gins 
7 


She he = 
{ | 


fens } -- 


; 
t 


Appendix F 


FLUID PROPERTIES 


F.1 Steam Properties 


Steam properties are one of the more essential parts of any thermal calculation in 
reservoir engineering. Interpolating data from the steam table needs the entire steam table to 
be included as an input file. Interpolation polynomials have been widely used to avoid this. 
Tortike and Farouq Ali (1989) presented a set of polynomials which were continuous and 
covered almost the full saturation envelope. This set of polynomials was used for the 
present work. 


F.1.1 Steam Condensate 


Viscosity 


27.1038 _ 23525.5 | 1.01425x10' _ 2.17342 x10" | 1.86935 x10" 


[ge = ~.0123274 + ep ete = z 


2791S STS 645:K (21) 


Density 


(,. = 3786.31 — 37.2487T + 0.196246T* — 5.04708 x 107 Fe 4629368 10a 
— 3.0848 x10°'°T? 


273.15 <T <640K (F.2) 
F.1.2 Steam Vapour 


Viscosity 


LL,, = -5.46807 x 10 + 6.89490 x 10°T — 3.39999 x10°T* + 8.29842 x Onde 
~ 9.97060 x10“ T* +.4.71914 x 10°77” 


273.15 <T< 645K (F33) 
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Density 


In p,, = —93.7072 + 0.833941T — 0.003208T? + 6.57652 x10°T? — 6.92747 x 10°T* 
a2 27 203.109 1 


Zi oA) = 1040s (F.4) 


F.1.3 Latent Heat of Vaporization 


L, = (7184500 + 11048.67 —88.40507? + 0.162561T? — 1.21377 x 10“T*)2 
273.15<T<645K (F.5) 
F.1.4 Saturation Temperature 
T, = 280.034 + 14.0856 In p + 1.38075(In p) —0.101806(In p) + 0.019017(In p)* 
611 kPa<p< 22.12MPa (F.6) 


F.1.5 Saturation Pressure 


Da (-175.776 9097270 1305 bao 0278 MAO Ee 


2973726 X10°T* +1.13816 x ee ye 


80 <T $647.3 K (ET) 


All the units in the above equations are in derived SI metric. Viscosity in Pa.S., 


pressure in kPa, density in kg/m3, latent heat of vaporization in kJ/kg and temperature in 


K. 


F.2 Oil Properties 


Viscosity 


Svrcek and Mehrotra (1988) used an ASTM viscosity-temperature chart D341 for 
seven different high viscosity petroleum samples to come up with a one parameter 
correlation. Although this is an adequate way to estimate viscosity of different bitumen, 
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their temperature range was limited to a maximum value of 130°C only. A good 


approximate method was suggested by Faroug Ali (1982), in terms of dynamic viscosity, 
in SI unit 


b/(T+273.15) 


H, = ae (F.8) 


If two points are known a and b can be easily found out. For the present study, two 
sets of a and b were found from the given values, 


a = 6.54053 E-15, for S33 oels 

b = 11380.31 (F.9) 
and 

a= 1.77632 E-7, for jeaew Eee « 

b = 4993.018 (F.10) 
Density 


A simple correlation was used to calculate the oil density using a thermal expansion 
factor, 


P. =P, x{1-B,(T -T.)}- (F.11) 
The subscript R represents a reference point (standard condition). 
Specific Heat 


Smith-Magowan et al. (1982) developed a correlation for the specific heat of bitumen 
from several samples of low, medium and high grade Athabasca oil sands. These were 
tested over a temperature range of 50 to 300°C, along with specific heats of components. 
Their equation gave an acceptable estimate of oil sands specific heat in J/°C g, 


c, = 1.557 + 5.219 x10°T - 8.686 x 10°T”. (F.12) 
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